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(ABSTRACT) 

A method has been developed which calculates two-dimensional, transonic, viscous flow in ducts. 
The finite volume, time marching formulation is used to obtain steady flow solutions of the 
Reynolds-averaged form of the Navier Stokes equations. The entire calculation is performed in the 
physical domain. The method is currently limited to the calculation of attached flows. 

The features of the current method can be summarized as follows. Control volumes are chosen so 
that smoothing of flow properties, typically required for stability, is not needed. Different time steps 
are used in the different governing equations to improve the convergence speed of the viscous cal- 
culations. A new pressure interpolation scheme is introduced which improves the shock capturing 
ability of the method. A multi-volume method for pressure changes in the boundary layer allows 
calculations which use very long and thin control volumes ( length/height = 1000). A special 
discretization technique is also used to stabilize these calculations which use long and thin control 
volumes. A special formulation of the energy equation is used to provide improved transient be- 
havior of solutions which use the full energy equation. 

The method is then compared with a wide variety of test cases. The freest ream Mach numbers 
range from 0.075 to 2.8 in the calculations. Transonic viscous flow in a converging diverging nozzle 
is calculated with the method; the Mach number upstream of the shock is approximately 1.25. The 
agreement between the calculated and measured shock strength and total pressure losses is good. 
Essentially incompressible turbulent boundary layer flow in an adverse pressure gradient is calcu- 
lated and the computed distribution of mean velocity and shear stress are in good agreement with 



the measurements. At the other end of the Mach number range, a flat plate turbulent boundary 
layer with a freestream Mach number of 2.8 is calculated using the full energy equation; the com- 
puted total temperature distribution and recovery factor agree well with the measurements when a 
variable Prandtl number is used through the boundary layer. 
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1.0 INTRODUCTION AND LITERATURE 
REVIEW 


Computational Fluid Dynamics ( CFD ) has become a powerful tool in the past 20 years in pre- 
dicting the fluid mechanics and heat transfer characteristics of certain flow fields. The advent of the 
modem, high speed computer has been the major driving force for the evolution and success of 
computational fluid mechanics. Because of the rapid decrease in computer costs as well as the 
equally rapid increase in computer speed, the solution of many complex fluid mechanics problems 
is possible on the digital computer. Computational fluid dynamics is used extensively in the aero- 
space and turbomachinery industries to predict flow over wings, around fuselages, through engine 
inlets, as well as the flow through components within the gas turbine engine. 

Computational fluid dynamics solves a fluid mechanics problem by first discretizing the flow do- 
main using a grid network of some kind. The governing partial differential equations are approxi- 
mated and an approximate form of the governing equations is used to solve for the discrete values 
of flow properties at node points of the grid network. There are a number of ways of approximating 
the governing equations (continuity, momentum, and energy) and the most common are the finite 
difference method, the finite element method, and the finite volume method. The net result of all 
these methods is that the non-linear partial differential equations ( governing equations ) are re- 
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placed by many simultaneous algebraic equations. These algebraic equations can then be solved 
for the unknowns, which are the flow properties at the node points within the grid. This discrete 
solution is an approximation of the continuous solution provided by nature. 


Experimental and analytical methods are also used in the analysis of fluid mechanics problems. 
Before the advent of the modem computer and modem computational algorithms, experimental 
methods were used extensively and were the main tool by which fluid systems were studied. The 
advantage of the experimental approach was and still is that it is capable of being the most realistic 
model of the true flow situation. However, the experimental approach does have a few disadvan- 
tages, namely ( 1) equipment is required to run the experiment (2) problems may arise in scaling the 
model to the actual flow, (3) measurement difficulties may be encountered ,and (4) the cost of the 
experiment may be excessive. 

Analytical methods have the restrictions that they are typically limited to simple geometries and to 
linear problems. They are usually restricted to linear problems because the mathematics involved 
becomes too unwieldy or the solution technique is not known. The exact solutions from analytical 
methods can be used to measure the accuracy of various numerical methods. However, because 
of the above restrictions on the geometry and linearity, analytical methods are of little practical use 
in solving today's complex engineering problems. 

When compared to the previously described experimental and analytical methods, numerical 
methods have the following advantages and disadvantages. Numerical methods have no restriction 
on the linearity of the governing equations. Because of this, complicated physical situations can 
be modeled. However, because of the discrete nature of the formulation, the solution accuracy is 
limited because of truncation error. This truncation error is caused by the coarseness of a grid and 
the type of discretization used in a particular algorithm. If the flow field is discretized using a finer 
and finer grid , the truncation error will typically be reduced. However, as the grid becomes finer, 
the cost of the computations may become excessive and the computer's storage capacity may be 
exceeded. Even though modem computational algorithms can be very powerful tools in the anal- 
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ysis of fluid mechanics problems, there are practical limitations to the problems which can be 
solved. 

The numerical techniques used to solve for the flow properties within a flow field can be classified 
into three main groups: they are (1) finite difference methods (2) finite element methods and (3) 
finite volume methods. The dominant method presently used in fluid mechanics is the finite dif- 
ference method. In the finite difference method, the continuous flow domain is 'discretized'' so that 
the dependent variables (flow properties) are considered to exist only at discrete points or consid- 
ered to the vary between grid points in a simple algebraic fashion (for example linearly). The gov- 
erning equations in differential form are approximated by replacing the derivatives in these 
equations with their approximate algebraic representations. The resulting algebraic equations are 
then solved for the unknown flow properties using an appropriate algorithm. 

The finite element method has been used extensively in the area of solid mechanics but it has not 
enjoyed as much success and popularity in the fields of aerospace and turbomachinery fluid me- 
chanics. In the finite element method, the flow domain is discretized into finite elements which are 
usually triangular in shape. The unknown properties are represented within each element by an 
interpolating polynomial which is continuous and derivatives of the properties are also continuous 
to a specified order within the element. As is typical of all finite element methods, information is 
obtained on an element by element basis and then this information is assembled into a global rep- 
resentation of the problem. To solve the governing equations, finite element solution techniques 
have used both the Galerkin method formulation and the variational formulation. 

A popular way of solving fluid mechanics problems in turbomachinery is the finite volume method. 
As with the finite difference and finite element methods, there are a number of different approaches 
which can be classified as finite volume methods. The identifying features of the finite-volume 
method are that the flow domain is descretized into finite control volumes and the governing 
equations in integral form* are used to solve for the unknown properties at grid points within this 
control volume network. When the equations of motion are written in integral form, the fluxes of 
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all quantities are identically conserved once the steady state is reached. If the finite difference for- 
mulation is used, all quantities are not necessarily conserved. Another advantage of the finite vol- 
ume method is its simplicity. Since the solution is obtained in a physical grid, no coordinate 
transformations are necessary and therefore the method is more easily understood. 

The finite volume method has been used extensively to solve the Euler equations for transonic flow 
including flow at high Mach numbers. Potential flow methods have been used in the past for cal- 
culations of inviscid transonic flow at Mach numbers close to one. However, potential flow 
methods assume that the flow is isentropic and therefore their solutions cannot satisfy the 
Rankine-Hugoniot shock relations. Therefore, when the Mach numbers in the flow field become 
greater than 1.2 to 1.3, the entropy increase across a shock wave becomes important enough that 
neglecting its effect would degrade the solution accuracy. For flows in which the Mach number 
exceeds 1.2 to 1.3, the Euler equations or Navier Stokes should be used so that the calculations 
correctly predict the shock strength and shock position within the flow field. However rather than 
solving just one second-order partial differential equation as is done for the two dimensional po- 
tential methods, the Euler equation formulation must solve the continuity equation, x and y mo- 
mentum equations, and the energy equation simultaneously. The finite volume formulation is 
preferred over the finite difference method for solving the Euler equations because of its property 
of conserving mass and momentum. 

When deciding upon a solution procedure which will calculate the properties at the node points in 
a transonic flow field, it should be noted that the steady form of the Euler equations is elliptic in 
nature when the flow is subsonic and hyperbolic in nature when the flow is supersonic. A different 
algorithm is required in the subsonic and supersonic regions of the flow field. The unsteady form 
of the Euler equations is always hyperbolic in nature; therefore the unsteady Euler equations may 
be solved with the same algorithm for both subsonic and supersonic regions of the flow. Shock 
waves in the flow field evolve as part of the solution. The solution to the steady problem is found 
by marching the unsteady solution in time until it reaches a steady state. This is the time marching 
method. 
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In the area of external aerodynamics, the principal development of the time marching, finite volume 
method has been by A. Jameson (1). In Jameson's method, as with all finite volume methods, the 
differential form of the governing equations is converted to the integral form using some form of 
Gauss' theorem. The unsteady form of the governing equations is used in obtaining a steady state 
solution by taking the asymptote of the unsteady solution. Jameson has solved the Euler equations 
for transonic flow over airfoils. The unsteady solution is advanced in time using an explicit 4-step 
Runge-Kutta method. He discretizes the flow domain into finite quadrilateral cells ( see Fig. 1.1 ) 
and places node points at the center of these cells. The properties at cell faces are determined by 
taking the mean of the properties at adjacent cell centers. This discretization results in a simple 
central difference scheme which must be augmented by the addition of dissipative terms which have 
a magnitude determined by the local flow gradients. These dissipative terms are designed to sup- 
press the tendency for odd and even point oscillations, and to limit the generation of wiggles and 
overshoots near shock waves. Currently Jameson has limited the use of his finite volume method 
to inviscid flow. 

In internal aerodynamics, McDonald ( 2 ) was the first investigator to use the time marching, 
finite-volume method. McDonald's method was limited to two-dimensions because of his choice 
of control volumes. Denton ( 3 ) extended McDonald's finite-volume method to three dimensions; 
the shape of the control volume was simplified and the solution procedure for the governing 
equations was modified to improve its accuracy. 

Versions of Denton's method have been used in inviscid viscous interaction programs for 
turbomachinery calculations. Singh ( 4,5 ) used Denton's method to calculate the inviscid region 
of the flow. The boundary layer was calculated using an extended Pohlhausen approach for the 
laminar flow. The turbulent region of the boundary layer was calculated using a method similar 
to that of Green ( 6 ) ( lag-entrainment ). Calvert (7,8, and 9) also used Denton's inviscid method 
for turbomachinery calculations. The laminar boundary layer was calculated using a momentum 
integral approach utilizing Thwaites' parameter. The turbulent boundary layer was calculated using 
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the method of Green et al. ( 10 ) which was further modified by East et al. ( 1 1 ) to handle regions 
of separation. The separation region was handled using an inverse boundary layer type calculation. 

Currently the principal workers using the finite volume method for internal flow field calculations 
are Denton(12,13,14), Van Hove(15), Holmes and Tong(16), Dawes(17) and Subramanian(18); all 
use the time marching approach. The majority of these finite volume methods have solved the 
Euler equations for inviscid transonic flow in turbomachinery blade passages [12,13,14,15,16 and 
18). Denton (13) has simulated the effect of the boundary layer using an inviscid- viscous interaction 
method. Fluid is injected through the surface of the blades to simulate the blockage effect of the 
boundary layer. An integral method is used to calculate the boundary layer. Dawes (17) has for- 
mulated a method which solves the Navier Stokes equations in a transonic compressor cascade. 
Each of these internal flow methods will now be discussed in more detail. 

The method of Holmes and Tong (16) is based upon the scheme used by Jameson. Modifications 
to the boundary conditions and the grid have been made to calculate flow in turbomachinery blade 
rows. The code is capable of calculating three dimensional inviscid flow in turbomachinery blade 
rows. The same fourth-order Runge Kutta scheme is used to advance the solution in time and the 
node points are again cell centered. Holmes and Tong comment that because cell centered node 
points are used, the accuracy of the scheme drops 'precipitously' when a non-uniform grid is used. 
A non-uniform grid results in control volumes of varying shapes and sizes. For inviscid calcu- 
lations, a flow field can be discretized such that the grid is fairly uniform. For viscous calculations, 
however, highly stretched grids are required and the accuracy of the scheme could become suspect. 
The method of Subramanian(18) also uses Jameson's finite-volume algorithm for inviscid 
turbomachinery flow calculations. 

Van Hove(15) has presented a finite-volume method for solving the Euler equations in a cylindrical 
coordinate system; this method is capable of calculating three dimensional flow in axial 
turbomachinery blade rows. Van Hove describes his method as a fully explicit, time marching, 
corrected viscosity, finite volume method. Artificial dissipation is needed to stabilize the transient 
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solution, but if the artificial dissipation were retained at the steady state the solution accuracy would 
be degraded. The corrected viscosity scheme is the means by which the artificial dissipation is re- 
moved as the steady state solution is approached. Typically some viscosity is retained to allow the 
capturing of shock waves without overshoots or undershoots. 

Of the finite volume methods used for turbomachinery calculations, the method of Dawes (17) is 
the only one to solve the Navier-Stokes equations. The grid points in the computation are cell 
centered like those seen in Fig. 1.1. The properties at cell faces are determined by assuming that 
the properties vary linearly between cell centers. This linear interpolation is contrasted with the 
method used by Jameson where properties at cell faces are the average of the cell centered proper- 
ties. The interpolation is needed for these viscous calculations because of the highly stretched grids 
used. A form of upwind differencing is used in the transient solution to enhance stability and is 
called a defect operator. This defect operator is added to the spatial flux operator to make the 
transient solution more robust and hence improve the stability of the method. The spatial flux 
operator is the actual flux imbalance for a given time level. To improve the accuracy of the sol- 
ution, this defect operator is reduced as the solution approaches a steady state. However, the 
convergence rate of the method is degraded when the centered differencing is retained near the 
steady state (see Fig. 1.2). The residual seen in Fig. 1.2 decreases rapidly when the defect operator 
is used; however, when the defect operator becomes small the residual appears to approach a steady 
state. This behavior is typical of the cell centered nodal structure when the artificial dissipation is 
removed. In addition, separate residual steady state artificial smoothing is needed. This residual 
smoothing is needed to capture shock waves and to prevent odd-even point oscillations. The 
method is implicit; however, the local time step cannot be unlimited and is typically set equal to 
10 times the explicit CFL condition. 
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Fig. 1.2 Convergence History for Dawes' Calculations 
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DENTON'S METHOD 


The method of Denton(12,13, and 14) has been widely used in the turbomachinery area. The new 
method to be discussed later in the body of this disseratation most closely resembles Denton's 
method. For this reason, his method will be discussed in more detail than the other finite volume 
methods. 

The control volumes are quadrilaterals cells with the node points located at the comers of the 
control volumes(see Fig. 1.3). This contrasts with the cell centered node arrangement used by most 
of the other workers. Properties for control volume faces are determined by taking the averages 
of the properties at the node points corresponding to the ends of a control volume face. It should 
be noted that Denton's method is not restricted to two dimensions but the following description 
of his method is based upon the assumption that we are dealing with a two dimensional coordinate 
system. 

In the previously described finite volume methods, the properties at node points are updated si- 
multaneously over one time step either explicitly or implicitly. Denton updates the properties se- 
quentially in the order of density, total energy per unit volume, pressure, (pu), and (pv) using the 
continuity equation, energy equation, equation of state, x-momentum equation, and y-momentum 
equation respectively. Properties are used from the previous time level to evaluate all terms in the 
governing equations except the updated pressure which is used in the momentum equations as soon 
as it is available. 

If the true pressure at a node point were used in the momentum equations , the solution procedure 
would be unstable. For this reason, Denton uses an effective pressure in the momentum equations. 
This effective pressure is interpolated from stirrounding grid points and the interpolation function 
does not include the actual grid point in its formulation. As the solution precedes towards a steady 
state, the effective pressure is adjusted towards the actual pressure. This adjustment is needed to 
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improve the accuracy of the final solution. However even at steady state, the effective pressure may 
not equal the actual pressure. This correction towards the actual pressure at steady state is another 
form of the defect operator used by Dawes. 

Denton's method is explicit in nature and he uses the CFL condition to determine allowable time 
steps. The CFL condition has been the criterion used in explicit methods to determine a maximum 
stable time step. Smoothing of flow properties is needed in Denton's method to stabilize the sol- 
ution. Typically this smoothing is only needed in the blade to blade direction. 

Since Denton's method does resemble the method currently under development, further details of 
Denton's method will be introduced in the body of the report. Denton's scheme will then be 
compared and contrasted with the current scheme. 


PRESENT CONTRIBUTION 


The scope of the present work is to extend a finite-volume method like that of Denton's to be able 
to calculate laminar or turbulent flow in ducts. The method will have the capability to calculate 
subsonic as well as transonic flow. In the process of developing the necessary ideas to extend 
Denton's inviscid method to viscous flow, the current work has also resulted in improved expla- 
nations for the inviscid methodology used by Denton. Some of the * art' used in Denton's 
method has been explained and explored. 

Some of the specific developments to be presented in the body of the report will be summarized 
here. Control volumes have been introduced which remove the necessity to smooth flow properties 
in the transverse direction. This development has significance in the reduction of numerical mixing 
and in the accurate modelling of viscous flows. A new pressure interpolation scheme has been de- 
veloped which results in solutions for transonic flow which have vastly improved shock wave defi- 
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nition. In contrast to all previous time marching methods where the same time step is used in all 
governing equations, different time steps are used in the continuity and momentum equations. 
These time steps result in improved convergence rates for viscous flow calculations. 

A method of discretization, called transverse upwind differencing, allows calculations to be made 
with highly stretched control volumes which are used in the viscous boundary layers near walls 
(aspect ratio [length/height] greater than 1000 ). A multi-volume approach for pressure changes in 
the boundary layer contributes to the stability of calculations with highly stretched control volumes 
and variable time steps. A new updating procedure is used which updates the pressure directly from 
the continuity equation. All previous finite-volume methods have updated the density through the 
continuity equation. 

A number of test cases are used to evaluate the accuracy of the new method. Test cases are also 
used to illustrate the effects of certain features of the method. For the moment, the test cases will 
just be listed. The details and purpose of each test case will be discussed later. Tes f cases 1 through 
7 use the assumption of constant total temperature as the energy equation. Test cases 8 through 
1 1 use the full energy equation. The test cases are 

1. Laminar boundary layer flow with a zero-pressure gradient. 

2. Inviscid transonic flow in Sajben's diffuser (23). 

3. Inviscid flow in a straight duct with inlet step profile in total pressure. 

4. One-dimensional nozzle calculations to demonstrate the new pressure interpolation scheme. 

5. Laminar flow in a converging duct with the boundary layer on one wall only. 

6. Turbulent boundary layer flow in an adverse pressure gradient approaching separation. 
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7. Attached boundary layer in transonic flow in Sajben's diffuser. 

8. Turbulent boundary layer in an adverse pressure gradient with an inlet freestream Mach 
number of 0.55. 

9. Rat plate turbulent boundary layer with a freestream Mach number of 0.95. 

10. Sajben's diffuser calculations including the energy equation. 

1 1. Flat plate turbulent boundary layer with a freestream Mach number of 2.8. 

Some of these cases will allow comparison with experimental or theoretical results ( cases 
1,2,4,6,7,10 and 11), some can be compared with other computational results ( cases 2,4, 5,6,7, and 
10 ) and some can only be judged qualitatively ( cases 3,8 and 9 ). 
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2.0 ANALYSIS SECTION 


2.1 GOVERNING EQUATIONS 


The unsteady form of the continuity equation, the x-momentum equation, and the y-momentum 
equation, in integral form, are used to obtain a steady-state solution for flow through 2-dimensional 
ducts by taking the limit of the unsteady solution as it approaches a steady value. The ideal gas 
equation of state and the assumption of constant total temperature complete the governing 
equations needed to solve for the unknown variables, p, u, v, P, and T. 

The unsteady continuity equation in differential form is 

-^- + V'pu = 0 . (2.1.1) 

at 

To transform the governing equations from differential form to integral form Gauss' theorem is 
used. Gauss' theorem says that 




( 2 . 1 . 2 ) 
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where ^ is an arbitrary vector and A is an outward normal area vector. When Gauss' theorem is 
applied to the continuity equation ( Eq. 2.1.1 ), we get 


+ V- p M ) dVol = dVol + jjpy • dA = 0 . (2.1 .3) 

For a finite control volume where we can assign one value of density to the control volume, and 
for a finite time step ,5 1, we get 


-|B-5Ko/= - Jjpy<U . (2.1.4) 

So for one time step in our calculation, we get a change in density ,8p, 

= (2.1.5) 

where the integral is evaluated explicitly at the current time step, n. The term in the brackets is the 
continuity error at time, n, therefore Eq. 2.1.5 can be rewritten as 

5p = (continuity errorW^rr * (2.1.6) 

oVol 

Previous finite volume methods , like Denton's (13), update the density after each time step using 
Eq. 2.1.6. When the density is updated in this way, the method will be referred to as a density 
update method. The current method , however, updates the pressure directly from the continuity 
equation for reasons which will be discussed later. When the pressure is updated directly from the 
continuity error, we will refer to the method as a pressure update method. In arriving at an ex- 
pression which relates the pressure change directly to the continuity error, we will assume that 
changes in temperature are small in comparison to other changes for one time step. With the above 
assumption about temperature changes, we can relate changes in pressure to changes in density 
through the ideal gas equation of state, in other words, 

5 P = RT5p. (2.1.7) 


2.0 ANALYSIS SECTION 


16 



If at this point the alternative assumption was made that the changes in pressure were isentropic, 
then 

5P = yRTdp . (2.1.8) 

If Eq. 2.1.6 and Eq. 2.1.7 are combined we get, 

5 P = { continuity error )RT- —~- . (2.1.9) 

oVol 

In other finite volume methods, like Denton's, the conservative form of the momentum equations 
is used to update the velocities for a control volume. The unsteady conservative form of the mo- 
mentum equations is 

dt~ + — ^ * P&y + V'\iVu+ V‘\iVy, r (2.1.10) 

where 8 , y is the Kronecker delta. When Gauss' theorem is applied to Eq. 2.1.10, we get 
(pw)” + 1 “ (pm)" = 5(p«) = 

( - Jfp UU'dA ~ tfPS,j-dA + Jj(n V U + \iW) ’ . (2.1.11) 

Another way of writing Eq. 2. 1 . 1 1 is 

8(pu) = ( momentum error )-^-r . (2.1.12) 

oVol 

The velocity at the new time step , n + 1 , is found by 

n + 1 _ (Pm)" + ‘ 
p" +1 

For the method introduced in the current work, a non-conservative form of the unsteady momen- 
tum equations is used. The non-conservative form is used because it allows the current method to 
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use different time steps for the continuity and momentum equations. The idea of different time 
steps will be discussed in further detail in section 2.5 ( TIME STEPS ). The* differences between 
the non-conservative form and the conservative form of the unsteady momentum equations are 
associated with the unsteady and convective terms. Specifically, we let 

+ V'puu = p^- + pa* Vu . (2.1.13) 

The right hand side of Eq. 2.1.13 can be rewritten as 


du 

IT + pa ' 


rj _ dU 

r “ =p lT 


+ V • p uu — n( V • pn) 


(2.1.14) 


In the current method, the right hand side of Eq. 2.1.14 is used in the momentum equation. 


Three observations can be made about the right hand side of Eq. 2.1.14. First, the velocity is up- 
dated directly from the momentum equation, rather than being updated indirectly through (pi/ ). 
Secondly, the second term in Eq. 2.1.14 is the conservative form of the convective terms used in 
Eq. 2.1.10. Finally, the third term can be recognized as the continuity error contribution to the 
momentum error. This term becomes zero at a steady state because the continuity equation at a 
steady state becomes 


F • pu = 0 . (2.1.15) 

When the new unsteady and convective terms are combined with the pressure and viscous terms 
and Gauss' theorem is applied, the momentum equations in integral form are 

(u) n+l - (U) n = 5(«) = 

l ~ ^PUU'd/i + Ufyu-d£ - JJP8,y*c/4 + JJXpPy + ^Vu r )-dA]r^~ ■ (2.1.16) 
Another way of writing Eq. 2.1.16 is 
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8u — ( momentum error ) — i*/ — . 

poKo/ 


(2.1.17) 


The ideal gas equation of state 


P = pRT 


(2.1.18) 


is used in other finite volume methods to update the pressure. In the current method the ideal gas 
equation of state is used to update the density since the pressure has already been updated through 
the continuity equation. Typically, the energy equation is simplified to the assumption of constant 
total temperature 

(u 2 + v 2 ) 

T 0 = T+ v L = constant . (2.1.19) 

2C P 


A more detailed explanation of the energy equation is included in a separate section ( ENERGY 
EQUATION Section 2.10). 


To maintain stability, the properties must also be updated in the proper sequence. In Denton's 
method, the sequence was 


1. Update the density from the continuity equation. 


2. Update the pressure from the equation of state. 


3. Update the (p«) and (pv) from the momentum equations using the new pressure and old ve- 
locities and old density. 


4. Update the temperature from constant total temperature. 


In the current method, the sequence is 
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1. Update the pressure from the continuity equation. 


2. Update the velocities from the momentum equation, using the new pressure and old velocities 
and old density. 


3. Update the density from the equation of state. 


4. Update the temperature from constant total temperature. 


2.2 CONTROL VOLUMES 


A new control volume has been introduced for this method. To eliminate the need for smoothing 
of flow properties, there must be as many control volumes across the duct as there are nodes where 
these variables are calculated. We need as many equations as unknowns. The control volumes also 
need to be located so that errors in continuity and momentum can correctly influence the changes 
in pressure or density and velocity without smoothing. The current control volume accomplishes 
this and is shown in Fig. 2.2.1. There are no nodes located along the wall. The nodes are located 
in the middle of the upstream and downstream faces of the control volumes. When calculating the 
flux through a streamwise face of an element, the values of the flow properties at the node on that 
face are used. This is equivalent to using step profiles on streamwise faces. When calculating the 
flux through a cross-stream face, first the properties are calculated at the comers of the element 
using linear interpolation, then the flux is calculated using the average of the flow properties at the 
ends of that face. 

Linear profiles on the streamwise face were also investigated. The results were essentially the same 
for linear and step profiles. Since step profiles are simpler they are used. 
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The control volumes used by Denton (13) look like those shown in Fig. 2.2.2. Fluxes of mass and 
momentum through each face are found by using averages of the flow properties stored at the ends 
of each face. However, since the number of nodes ( unknowns ) is greater than the number of 
control volumes ( equations ), smoothing of flow properties is needed in the crossflow direction to 
remove the dependence of the final solution on the initial guess. 

If the current method were to be extended to three dimensions the following additions would need 
to be made. An example of a typical three dimensional control volume is shown in Fig. 2.2.3. 
The location of control volume boundaries are specified in the input data and the control volume 
surfaces are constructed from this information. Once the control volume boundaries are known 
then the grid points are placed in the middle of the upstream and downstream faces of the control 
volume. The fluxes though the transverse faces of the control volume needed for the continuity 
and momentum balances are determined from interpolated properties using the nodes adjacent to 
the face. Fig. 2.2.4 shows two adjacent control volumes of different sizes. The procedure for cal- 
culating the properties to be used in calculating the fluxes on the common boundary ( face I ) can 
be shown in the following way. For face I, any property, X, is determined from the average of the 
property at points A and B, where the values of the properties X A and X B are determined by linearly 
interpolating between the values of the property at nodes 1 and 2, and between the values of the 
property at nodes 3 and 4, respectively. 

Assuming that face II corresponds to a solid boundary, the values of a property at points C and 
D are determined by linear extrapolation using the values of the property at nodes 1 and 2, and 3 
and 4; respectively. For the present calculations, only the pressure needs to be calculated at the 
solid boundaries since the fluxes of mass are set equal to zero through these solid boundaries. 
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Fig. 2.2.4 New Control Volumes in Non-Uniform Grid 



2.3 DISTRIBUTION OF PROPERTIES 


The properties at node points are changed in the flow field after each time step because the conti- 
nuity and momentum equations are not satisfied for a given control volume. The amount that 
properties are changed at nodes depends upon the extent to which continuity and momentum are 
not satisfied, the volume of the control volume, and the time step. A decision must also be made 
about which node, either upstream or downstream ( I or I + 1 in Fig. 2.3. 1 ), these changes should 
be allocated to. The criterion to be used in determining where changes in properties should be sent 
is that these distributions result in reduced errors in continuity and momentum. Let us start with 
the momentum equations. 

In the current method, the momentum equation specifies a change in velocity for the control vol- 
ume based upon the momentum error ( Eq. 2.1.17 ), 

Su = ( momentum error ) — - . (2.3.1) 

po Vol 

This change in velocity is assigned to the downstream node of the control volume. This change 
will result in a stable calculation procedure and can be explained using a simple example. 

Refering to Fig. 2.3.2, if we assume that the only contribution to the momentum imbalance is 
caused by the convective terms through the streamwise faces, the momentum error for that control 
volume is 


momentum error = (p jUjAfiuj — (p i+]U I+ \A l+ \)u I+ 1 . (2.3.2) 

Suppose that the net momentum error for the control volume is positive( in other words, 
(p iU,A t )u, > (p /+1 « /+1 /l /+1 )u /+1 ). Using Eq. 2.3.1, this error will result in an increase in velocity for 
the control volume and this increase will be assigned to the downstream node(I + 1). This increase 
in velocity at I + 1 will act through Eq. 2.3.2 to reduce the momentum error for the control volume. 
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Therefore sending the velocity change to the downstream node is stable. If the change in velocity 
(increase) was sent to the upstream node , the momentum error would increase and this distribution 
would therefore lead to an unstable procedure. 

When the conservative form of the momentum equation is used by Denton ( Eq. 2.1.10 ), the x- 
momentum error results in a change in (pu) for a control volume. For the same reasons as outlined 
above, the change in (p u) is also sent to the downstream node. 

In the current method, the pressure is updated directly from the continuity equation ( pressure 
update method ) rather than the density being updated ( density update method ) . The pressure 
change is assigned to the upstream node ( I ) ( see Fig. 2.3. 1 ) to satisfy stability considerations. 
The change in pressure for a control volume is related to the continuity error by 

8 P = ( continuity error )R T-zjy t . (2.3.3) 

oVol 

Refering to Fig. 2.3.3, if the error in continuity is due only to errors in streamwise velocity and 
density then the continuity error for that control volume is 

continuity error = p jUjAj - p /+1 u /+1 i4 /+1 . (2.3.4) 

If, for example, the continuity error is positive for the control volume, the pressure will increase for 
the control volume ( see Eq. 2.3.3 ) and that increase will be assigned to the upstream node. Let 
us assume for clarity that the momentum equation is initially balanced. This pressure increase will 
change the surface force on the face of the control volume coincident with node I and this force 
will unbalance the momentum equation and introduce a momentum error for the control volume 
of interest and also for the upstream control volume. This increase in pressure at node I will result 
in a positive momentum error for the control volume of interest and a negative momentum error 
for the upstream control volume. The velocity at node 1+ 1 will increase ( 5 u l+l ( + )) and the 
velocity at node I will decrease ( 5w 7 ( — )) because of these momentum errors. These changes in 
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-D Example With Continuity Imbalance ( Pressure Update Method 




velocity reduce the error in continuity ( see Eq. 2.3.4), therefore sending pressure changes upstream 
is a stable procedure. 

In contrast, if the pressure change was sent downstream, the momentum equation would cause the 
downstream velocity to decrease. This distribution causes the continuity error to increase and this 
decision would result in an unstable calculation method. 

In the current method, the density used in the continuity and momentum equations may be dif- 
ferent than that determined from the ideal gas equation of state. The reason for using an effective 
density rather than the true density is that the density change at a node will be of the same sign as 
the pressure change since bP oc 8p RT. These density changes may act to violate the stability cri- 
terion that the distribution in properties results in improved conservation of mass and momentum. 
Therefore an effective density which satisfies our stability criterion must be chosen. This topic will 
be discussed in more detail in section 2.4 (PRESSURE INTERPOLATION METHOD). 

In the finite volume method of Denton , if the continuity error for a control volume is not zero, 
the continuity equation specifies a change in density for the control volume (density update 
method) 


5p = (continuity error) ■ — r (2.3.6) 

o Vol 

and the change is typically sent to the downstream node ( 1+ 1 ) shown in Fig. 2.3.1. This dis- 
tribution of density results in a stable calculation procedure because of the following reasoning (refer 
to Fig. 2.3.4). If, for example, the continuity error is positive, in other words, too much mass flows 
into the control volume, then the continuity equation will specify an increase in density ; Eq. ( 2.3.6 ), 
If the density at the downstream node is increased (5p /+1 ( + )) then this increase in density 
will act through Eq. 2.3.4 to decrease the continuity error. So sending density changes to the 
downstream node will result in a stable method. In contrast, if the density change was sent up- 
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stream ( Sp f ( + ) ) then the continuity error would increase. This decision would result in tin 
unstable calculation method. 

When Denton uses the density update method and sends changes in density to the downstream 
node, the pressure used in the momentum equations is an effective pressure which may be different 
from the pressure determined from the ideal gas equation of state using the density and temperature 
at that node. This is because the density change at the downstream node will cause a pressure 
change at that node in the same sense as the density change. This pressure change, acting through 
the momentum equation, will cause the velocity to change at the downstream node and this velocity 
change wiJl contribute towards violating the stability criterion. For example, the positive continuity 
error discussed above would result in an increase in pressure at node 1+ 1. This pressure change 
would cause the momentum equation to have a negative momentum error. The downstream ve- 
locity would decrease because of this error and a decrease in velocity at node I + 1 (see Eq. 2.3.4 ) 
would increase the continuity error. It is because of this interaction between density changes and 
pressure changes that the density update method must use an effective pressure in the momentum 
equation which will satisfy our stability criterion. This topic will be covered in more detail in sec- 
tion 2.4 ( PRESSURE INTERPOLATION METHOD ) . 

In review, the current method uses the following sequence in updating properties over one time 
step. 

1. the pressure is updated through the continuity equation and the pressure change is sent to the 
upstream node. 

2. the u and v velocities are updated through the momentum equations using the new pressure 
and the velocity changes are sent to the downstream node. 

3. the effective density is updated through the ideal gas equation of state using an interpolated 
pressure. 
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4. the static temperature is updated from the assumption of constant total temperature. 

In the finite volume method of Denton, the following sequence is used to update properties over 
one time step. 

1. the density is updated through the continuity equation and the density change is sent to the 
downstream nodes. 

2. the pressure is updated using the ideal gas equation of state and the updated density. 

3. an effective pressure is calculated from an interpolation formula using the new pressures. 

4. (pu) and (pv) are updated through the momentum equations using the new effective pressure 
and changes in (pu) and (pv) are sent to the downstream nodes. 

5. the static temperature is updated from the assumption of constant total temperature. 

Both the sequence in which properties are changed and the nodes to which these changes are sent 
are very important in maintaining stability for these finite volume methods. 


2.4 PRESSURE INTERPOLA TION PROCEDURE 


As part of the updating procedure used by Denton (13), an effective pressure is used in the mo- 
mentum equations rather than the true thermodynamic pressure determined from the equation of 
state. This effective pressure is needed because if the true pressure is used in the momentum 
equations the solution may not converge. A simple explanation for this instability is that pressure 
changes caused by the continuity error will act through the momentum equation to cause the error 
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in continuity to increase rather than decrease. As discussed in the section 2.3 (DISTRIBUTION 
OF PROPERTIES) this will lead to an unstable calculation method. 

In the current method, the density used in the continuity and momentum equations is an effective 
density which may be different than the density obtained using the ideal gas equation of state. The 
use of an effective density is required because of stability considerations. If the actual density is used 
in the continuity and momentum equations, the solution procedure may be unstable as mentioned 
in section 2.3 (DISTRIBUTION OF PROPERTIES). 

Stability of Denton's Method: In order to achieve stability for Denton's scheme the pressure used 
in the momentum equations is not the true pressure but an effective pressure, P which is only 
an approximation to the true pressure. The use of the effective pressure is described by Denton (13) 
in the following manner. Stability is obtained by setting 

PeffJ = ^/+ 1 + CF, (2.4.1) 

where CF, is a correction factor which is an approximation to the difference in pressure between 
node points I and I + l(see Fig. 2.4. 1). The scheme is' unstable if CF, is found directly from the 
difference in the true pressures, even if the changes in CF, are highly damped in time. However, 
several stable methods of estimating CF, are possible and two of them are used by Denton. With 
reference to Fig. 2.4.1, a simple estimate of the correction is 

CF, = ~ (2.4.2) 

where a = 1 would give a centered difference 1st order accurate estimate of the correction. It 

should be noted that the accuracy of the solution is determined by the accuracy of the difference 

in effective pressures between adjacent grid points and hence the above estimate will give a second 
order accurate solution. 
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MERIDIONAL DISTANCE, m 


Fig. 2.4.1 Calculation of Effective Pressure 
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An alternative scheme is to base the correction on a parabolic interpolation between points 1+ 1, 
1 + 2, and 1-1, This leads to 

CFj = ~ P,+2) (2.4.3) 

which is second order accurate when a = 1. In order to provide sufficient numerical damping to 
capture shock waves cleanly it is found desirable to make the value of a slightly less than 1, say 0.8 
-0.9. An optional alternative is to make the value of a depend on the density gradient in such a 
way that a is automatically decreased in the region of a shock wave. The expression for a , based 
upon the density gradient is, 


a = (1 - . (2.4.4) 

This method provides additional stability in regions of large gradients in properties like those seen 
through a shock. The density gradient correction is applied only when there is an increase in den- 
sity. 

Stability of the Pressure Update Method: Currently the pressure update method is used and an 
effective density based upon an interpolated pressure can be used to stabilize the calculation pro- 
cedure. Using a three point interpolation scheme, like Eq. 2.4.3 above, one possible expression for 
calculating an effective density is 


P/+1 = I P, + — t P/ ~ 2) l x • (2.4.5) 

This effective density will result in a stable calculation procedure for all Mach numbers but will 
smear a shock wave out over several grid points. 

As mentioned previously, when the effective pressure or effective density used in the governing 
equations is relaxed to ideal gas, the solution procedure may become unstable. In the remainder 
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of section 2.4, this instability mechanism will be investigated, the above three point interpolation 
schemes will be shown to be stable, and a new Mach number dependent interpolation scheme will 
be introduced which gives much better shock capturing properties than the previously used two or 
three point interpolation schemes. A 1-D in viscid flow example will be used as a vehicle to develop 
these ideas. These ideas were developed by J.G. Moore in reference 19. 

1-D Flow Example: For the geometry shown in Fig. 2.4.2, we are seeking a 1-D steady flow sol- 
ution which satisfies the continuity equation 


V • pu = 0 


and the momentum equation 


V * PU U = - VP 


(2.4.6) 


(2.4.7) 


and which also satisfies the ideal gas equation of state and maintains constant total temperature 
throughout the flow field. 


Continuity: In discretized or integral form, the continuity equation for a converged 1-D solution 
becomes, 

P/+1 N/+1 ^/+i- p/U/i4/ — 0 (2.4.8) 

where the superscript f stands for the final converged values. 

For an intermediate solution the velocity will be u and the density will be p. The differences be- 
tween the current properties and the correct properties are 5 u and 5p. The continuity equation, 
Eq. 2.4.8, can be rewritten as 

(P/+1 + SP/+i)( u /+i + iMz+i “ (P/ *+* 5p/)(w/ + SufiAj = 0. (2.4.9) 

Expanding and rearranging, Eq. 2.4.9 becomes 
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P/+l^/+i^ u /+l — P lAfiui + «/+i./4/ + i5p| +1 UjA/dpj — 


{PjUjAj P/+i“/+i^/+ 1 ) + (5p/5u jAj 5p/ +1 5u /+1 /4 /+1 ). (2.4.10) 

The first two terms on the right hand side of Eq. 2.4.10 represent the continuity error and the last 
two terms are of the order & and will be negligible when the computation is nearly converged ( i.e. 
5p < < p and 5 u < < u ). 

Therefore Eq. 2.4.10 can be rewritten as 

P/+l^/+i 5u /+i ” P/^/H + u,+ 5p /+1 - UjAjSpj = m erwr + small. (2.4.11) 

In a time marching method, the left hand side of Eq. 2.4. 1 1 can also be used to evaluate the change 
in mass flow rate for one time step where the changes in properties in Eq. 2.4. 1 1 are for one time 
step. 

In the density update time marching calculation procedure (Denton), the density is updated from 
the continuity error. From the integral form of the unsteady continuity equation, the density on 
the downstream side of the control volume is changed according to 

8p I+x =m error4 -$L- . (2.4.12) 

The density change, 5p /+lJ affects the change in mass flow rate directly, but it also acts through the 
ideal gas equation of state to change the pressure. This pressure change then acts through the 
momentum equation to change the velocity ( 5 u /+l ). 

In the pressure update time marching calculation procedure ( Nicholson/Moore ), the pressure is 
updated directly from the continuity error. The pressure is changed at the upstream node according 
to 
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5/> / = terror J^ t c RT l Vol l • 


(2.4.13) 


This pressure change acts through the momentum equation to change the velocities, u /+l and u„ 
and it also acts through the ideal gas equation of state to change the density. We will now look 
at the discretized form of the momentum equation for our 1-D example. 

Momentum: The steady state momentum equation discretized over the control volume between 
points I and I + 1 is 


(P u A)j+\ «/+i (puA)[Uj- Ify+i^z+l PjAj P s (A t+ i Aj)\ (2.4.14) 

where P, is the pressure acting on the sides of the control volume. Traditionally, the pressure has 
been 


_ (/>/+! + Pj> 

2 


(2.4.15) 


Eq. 2.4.14 can therefore be rewritten as 


(P wy O/+i w /+i (p uA)jUf— (P { + j Pi)(Ai+ j + Aj)/2. (2.4.16) 


We may rewrite Eq. 2.4.16 as, 

" I /+i u ;+i “ ~ ~ ( p i+\ ~ Pi)Vol,/8x, (2.4.17) 

where m — p uA is the local mass flow rate, Vol, = 8x^A /+l + A,)/2 is the volume of the control 
volume and Sx, = x /+1 — x t is the grid spacing for the control volume, I. Eq. 2.4.17 may be re- 
written as 

( u /+l ~ «/)(«/+ 1 + "»/)/2 + (m /+ , - m,)(u I+l + u,)/ 2 = - (P /+1 - Pj)Vol,/5x, (2.4.18) 
or 


2.0 ANALYSIS SECTION 


41 



"tMi + 1 ~ «/) + um errorJ = - (P /+1 - P,) Vol,/5xj . 


(2.4.19) 


In the current method, the continuity error term , um^j , is omitted because of the non- 
conservative form of the momentum equation that is used ( Eq. 2.1.16) . Therefore, the change in 
velocity on the downstream side of the control volume is changed in direct proportion to the mo- 
mentum error. From the unsteady form of the momentum equation, the change in velocity at the 
downstream node is 

5 U/+, = I -(/*/+! + 5 />,+ , - Pi - 8 P,)Vol,lSx, - m(u l+l - u t m{p l+l Voi t ) (2.4.20) 
where the pressure change, 8 P, has been calculated from the continuity error. 

In the method used by Denton, the conservative form of the momentum equation is used to cal- 
culate the change in (pu) ;and that change is 

5(P«)/+| = u /+ ,Sp /+I + p /+l 8u /+ , = 

[ — (Py+l + 8P l+ j — Pj — 8P j)Volj/5xi — rn{u I+ j — uf) — Um error j\8tlVoi l . (2.4.21) 

By taking the mean velocity for the control volume , w, approximately equal to the velocity on the 
downstream side , u l+u we may subtract u times Eq. 2.4.12 from Eq. 2.4.21 to get the change in 
velocity, 8u, + 1 , for one time step. The resulting equation is 

5u /+i = l ~ ( p l+\ + 5/> /+i - Pi “ 8Pj)Vol,l8x } - m(u I+i - u,)\8tl(p, + i Vot,) (2.4.22) 

having assumed that the same time step is used in both density and velocity updates. This is the 
same result that is obtained using the Nicholson/Moore method ( Eq. 2.4.20 ). 

Let us assume that at the beginning of a time step that the momentum equation is balanced except 
for the pressure change introduced through the continuity error, in other words, 
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Mu I+ , -U,)= - (P /+1 - P,)Voljl5x,. 


(2.4.23) 


The momentum balances for both methods, Eq. 2.4.20 and Eq. 2.4.22, give us 

5 u l+ , = (8 P, - 5P l+ , )-- 5 ' ■ . (2.4.24) 

(P/+ i5*/) 

In general, in the density update ( Denton ) method, the time step is calculated from the CFL 
condition, which for 1 -dimensional flow is 


8/ £ — (2.4.25) 
(u + c) 

where c is the speed of sound. In the pressure update method the time step for momentum is ob- 
tained from the coefficient of « /+I in the steady flow momentum equation and is 

8f = -^- . (2.4.26) 

A more detailed description of the time steps is presented in section 2.5 ( TIME STEPS ) . We 
may combine these two equations by saying 


8f = — — — (2.4.27) 

(m + ec) 

where e = 1 for the density update method and e = 0 for the pressure update method. Eq. 2.4.24, 
which is valid for both methods becomes, 


. ... _ ( 5 />, -«>,+ ,) 

P,+ ,8 “' +l_ („+«) 


(2.4.28) 


The changes in density and velocity can now be substituted into Eq. 2.4. 1 1 to evaluate the change 
in mass flow rate for one time step. 
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Change in Continuity for One Time Step: The left hand side of Eq. 2.4. 1 1 may be used to evaluate 
the change in mass flow rate for one time step. Substituting Eq. 2.4.28 into Eq. 2.4. 1 1 to eliminate 
p8u yields 


A I+ 1 ~r 1 , . /+1 ^ A l “ 7~~ r ~ + u I+\ A r+\&Pl+\ ~ u I A fiPl ~ ™changej (2.4.29) 

(u + ec ),+ , (u + ec), 


If we now rearrange Eq. 2.4.29 to order the coefficients of the 8 P's and 8p'r , we get 

-A, 


(u + ec); 

I A/+j 
(u + ec);+] (u+ec); 

A, 


5/V 


/-l 


(u + ec ),+ , 


+ 77-T77, 1 5P l ~ u l A fiPi 

5P l+ i+ u /+1 i4 /+1 8p /+1 ] 


^ change ,/ • (2.4.30) 


For stability we require that the change in mass flow rate be of the same sign as the error in mass 
flow rate. Note that this stability requirement is a necessary condition for stability but it may not 
be a sufficient condition to ensure stability. 


Stability of the Density Update Method Using Ideal Gas: The density update method updates the 
density directly from the continuity error and then the pressure is updated through the ideal gas 
equation of state. For an intermediate solution where there is a continuity error only between 
nodes I and 1+ 1, Eq. 2.4.12 yields 


8p /+ i - m error , 


8/ 

Volj 


(2.4.31u) 


8p ; = 0 (2.4.318) 

8p;_i = 0 (2.4.31c) 


and from ideal gas, assuming temperature changes over the time step are negligible, 
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(2.4.32a) 


5/ 


& p i + 1 ~ 8p/+i/?T - nt error jRT-j^j- 


8 P, = 5p ,RT = 0 


— 5p — 0. 


Substituting these results into Eq. 2.4.30, we get 


r * ^/+ j RT ^ 

(u + c) W ^ + l / */ + 1 \ m errorJ y 0 [^ = m change,l 


Since for stability we require that fh,, re , and have the same sign, we must have 


I <^7 + “" 1 |^'> 0 - 


If we substitute c^/y for RT into Eq. 2.4.34, we get 


— c 


y(u + c) /+ , 


+ «/+ 1 >0 


or 


Y« /+ i(u + c) /+1 > c 


For a y of 1.4, Eq. 2.4.36 will be satisfied if 


u > 0.48c 


or 

M > 0.48 . 

Thus for low Mach number flows, this density update method is unstable. 


(2.4.326) 

(2.4.32c) 

(2.4.33) 

(2.4.34) 

(2.4.35) 

(2.4.36) 
(2.4.37a) 
(2.4.376) 
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Stability of the Pressure Update Method Using Ideal Gas: The pressure update method updates 
the pressure directly from the continuity error and then the density is updated through the ideal gas 
equation of state. Again, for an intermediate solution where there is a continuity error only between 
nodes I and 1+ 1, Eq. 2.4.13 yields 


SPl = menoi-focRT/Volf 

(2.4.38a) 

5P /+ , = 0 

(2.4.386) 

03 

II 

o 

(2.4.38c) 

and from the ideal gas equation of state assuming that temperature changes over one time step are 

negligible, we get 


x 8P, St e 

5p/ rt 

(2.4.39a) 

- o . 

(2.4.396) 

Substituting these results into Eq. 2.4.30, we get 


/A/+\ /l i tifAj , $t c RT 

( U /+ , + U, RJ ) m error,l y 0 ^ m change,I- 

(2.4.40) 

For m„ nr and to have the same sign, we will require that 


Aj + 1 + A r _ UjAj ^ 
“/+! u l RT 

(2.4.41) 


If we assume that the values of properties at node I are approximately the same as the values at 
node 1+ 1 and if we substitute c*/y for RT, Eq. 2.4.41 becomes 



(2.4.42) 
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or 


u < 



(2.4.43) 


For y = 1.4, for stability to be maintained 


u < 1.2c (2.4.44a) 

M < 1.2 . (2.4.446) 

Thus for high Mach number flows, this pressure update method is unstable. 

A Downwind Effective Pressure or an Upwind Effective Density Method: If an inconsistency in the 
pressure-density relation ( the ideal gas equation of state ) is introduced such that the pressure used 
in the momentum equation is offset by 1 grid point from the density used in the continuity 
equation, the equation of state may be written as 

Pi=p l+l RT (2.4.45) 

In a density update method, the pressure used in the momentum equation, an effective pressure , 
is evaluated using properties downwind of the actual node. Similarly, in a pressure update method, 
the density used in both continuity and momentum equations is evaluated using a pressure upwind 
of the actual node. For both the density and pressure update methods, the changes in density and 
pressure are 

ommi 


8p/-i = 0 


(2.4.46a) 


£ 01 

0P[ — n^error,l-\ p 


(2.4.466) 
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^P/+ 1 ™ err or, I 


bt 


VoL 


7+1 


(2.4.46c) 


ERESSUM 


bt 


bP/-i m error j-yRT y o i^ ^ 


(2.4.46d) 


bt 


M/ = terror, /RT-y^ 


(2.4.46e) 


bPi + 1 n i err orj+ \ PT- Yq! 


bt 


7+1 


(2.4.46/) 


Eq. 2.4.30 now becomes, 


_ . Ai . UjA. StRT v 
\u+ec), RT )K Volj. x } 


m errorJ - 1 


+ / dill + di + “Z+l-^Z+l y 8f#7\ 
(u + ec)i + 1 (« + ec){ RT Vol { 


™errorJ [ ^change ,/ * (2.4.47) 


-(■ 


+ 1 w StRT \ 

(u+ec) /+1 Vol t+ , 


m error,I+ 1 


From this equation we can see that the coefficient of m €rror I is always positive and so the downwind 
effective pressure method and the upwind effective density method both pass the simple stability 
criterion ( has the same sign as ) for all Mach numbers. It should also be noted that 

the coefficients of ffi trror%I _ 1 and rh mt%l + j are of the opposite sign to the coefficient of m errorJ and it 
is generally of a smaller magnitude. This further assures the stability of Eq. 2.4.45. 


While the pressure-density relationship of Eq. 2.4.45 is stable, testing has shown that it results in 
solutions with poor shock capturing. The calculated shock is spread over numerous grid points. 


2.0 ANALYSIS SECTION 


48 



Fig. 2.4.3 shows the calculated and theoretical pressure distributions for a l-D calculation using 
Eq. 2.4.45 with a nominal shock Mach number of 1.45. 


Stability of a 3-Point Interpolation for the Effective Pressure: One of the pressure-density relations 
which can be used in the density update method is a three point interpolation of the density used 
in the ideal gas equation of state to obtain an effective pressure. If we assume that the temperature 
is approximately uniform, we may write the interpolation formula as 

P, = (P/ + i " {Pl + 2 ~ Pl ~ t) ) RT . (2.4.48) 


The change in pressure for one iteration is therefore 


5P, = (5p /+1 


Sp /+2 

3 



(2.4.49) 


Substituting this change of pressure into Eq. 2.4.30 and neglecting variations of A,u, and c with x, 
we obtain, 


/ Ac 2 s/s;„ Sp /+t 8 p/- 2 N 

~ ( 77 — rr)(5p/ - — T — + — x — ) 
(y(u + c)) 3 3 


+ (- 


2Ac l 


(y(u + c)) 




= m. 


1 change ,/ 


(2.4.50) 


- (- 


Ac 


(Y (u + c)) 


■)(Sp /+2 ~ 5 p l +3 + ~r~) - uAS p /+1 


Collecting terms of 5p and substituting the Mach number , M , for u/c , we get 
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(2.4.51) 


+ Acl( 3y(A/ + 1)) Sp /+3 

- SAcHMM + 1)) 5p /+2 

+ Ac{7/(MM + 1)) + M) 8p /+1 

- /1c(4/(3y(M + 1)) + M) 8p, 

+ 2Ac/(3y(M + 1)) 8p,_ , 

- Ac/(3y(M + 1)) 8 P /_j 

= m change,l • 

We can see that the change in mass flow rate for the control volume between grid points I and I + 1 
is now dependent upon the the changes in density at 6 grid points. If there was only a mass flow 
rate error for control volume (IJ) the simple test for stability would be satisfied. Since the coeffi- 
cient of 8p /+1 is positive for all Mach numbers the change in mass flow rate would have the same 
sign as the error in mass flow rate. 

However, since the coefficients of 8p /+J and 8p,_! are also positive, a more sophisticated stability 
criterion is appropriate. The stability criterion that we will apply is : the center point coefficient 
must be greater than the sum of the other positive coefficients. 


£( centerpoint coefficient ) > Z( other positive coefficients ) 


(2.4.52) 


Applying this criterion to Eq. 2.4.51, stability requires that 


Ac(7/(2y(M + 1)) + M) > 3Ac/(3y(M + 1)) 


(2.4.53) 


This expression is always true, therefore Eq. 2.4.48 should be stable for all Mach numbers. The 
experience of Denton and other users of his code confirms this. 


Stability of the 3-Point Interpolation for the Effective Density: A similar analysis can be done for 
the pressure-update effective-density method using a three point interpolation of pressure to obtain 
the effective density 


2.0 ANALYSIS SECTION 


51 



P/+ 1 ~ [ p i + 


(2.4.54) 


( P l + 1 ^/-2)i 1 

3 RT ' 


For the change in density, we get 


8p/+ 1 = (5 P[ + 


3 'RT 


(2.4.55) 


and then substituting this change into Eq. 2.4.30, gives us 


+ (Alc)( - 1/A/ + yM/3) 5 P I+ , 

+ {A/c){2/M + 2yA//3) 5P, 

- {Alc)(l/M + yM) SPt-i 

- (A/c)(yMI3) SP /_2 

+ {Alc)(yM/3) 5P/_ 3 


m change,I • 


(2.4.56) 


The center point coefficent is the coefficient of 5 P, , since this is proportional to m„ rorJ in the 
pressure update method. In Eq. 2.4.56, the coefficient of 5 P, is positive and greater than the sum 
of the other positive coefficients; therefore Eq. 2.4.54 should be stable for all Mach numbers. 


Mach Number Dependent Interpolation Formula for Effective Density: From experience, it has 
been observed that when the Mach number is low, the pressure update method is stable with the 
ideal gas equation of state satisfied at each grid point. Since this is the correct pressure-density re- 
lation for ideal gases it should be used where feasible. In this section we will start with a generalized 
pressure interpolation equation for the effective density 

P/+i - \ p i + a o( p i + 1 - p i ) + a \ 2 + 02 3 '~RT (2A57) 

and seek Mach number limitations at and Oj using criterion ( Eq. 2.4.52). Comparing 
equations 2.4.45 and 2.4.57, the upwind effective density corresponds to q, = a x = Oj = 0, the three 
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point interpolation, Eq. 2.4.54, corresponds to q, = a, *= 0, a, = 1 , and ideal gas to 
q> = l.a, = Oj = 0. Substituting 

5p /+1 = |(1 - q^P, + (q, + ^- + -y-)8/> i+ j - ^P,-x ~ ^ZP t - (2.4.58) 

into Eq. 2.4.30 and rearranging in terms of the coefficients of each 8 P, q,, a u and a, yields 


(4)[( - MM + y M^ + (yM/2)a x + (yM/i)a 2 )5P I+ , 

+ (2/M + yM - 2yM% - (yM/2)a x - (yM/3)o 2 )8/> / 

(-\/M-yM + yMaQ - ( yM/2)a x )bP f _ ! 

( + (yM/2)a x - (yM/l^SP,^ 

) 

+ ( + (yA//3)a 2 )8P / _ 3 ] 

Let us first consider the case when a, = q, = 0 and find limiting values of q,. From Eq. 2.4.57, it 
is obvious that we should consider only values in the range 

0 £ q, £ 1 . (2.4.60) 

The coefficient of 5P, is positive when 

2/M + y M - 2yMoQ > 0 . (2.4.61) 

This gives a limit on q, which is a function of the Mach number and the limit is 

qj < 1 /(yM 2 ) + 1/2 . (2.4.62) 

But the coefficient of 5P /+1 is positive when 

- I/M + y Moq > 0 , or M 2 > \/(y q,) . (2.4.63) 

In this region , from Eq. 2.4.52, we require that 


> — m, 


'change,! 


(2.4.59) 
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2/A/ + yM — 2y M oq > — 1/A/ + yA/og 


(2.4.64) 


or 

q, < l/(yA/ 2 ) + 1/3 . (2.4.65) 

Valid values of q, based on these criteria are shown as a function of Mach number in Fig. 2.4.4. 

Let us next consider limiting values of a, when q, and dj are zero. From Eq. 2.4.59 the coefficient 
of 5 P, 

21 M + yM - yA/ 0,/2 > 0 (2.4.66) 

is positive for all Mach numbers in the range 

0 £ a, ^ 1. (2.4.67) 

The coefficient of SPj -2 is positive for all M and the coefficient of 8 P I+l is positive when 

- 1/A/ + yMa x /2 > 0 or M 2 > 2/(y a t ). (2.4.68) 

For M 2 < 2/(y a,) we then require the coefficient of 5P, to be greater than the coefficient of 8P,. 2 > 

2/ M + yM - yMa x /2 > yMa x /2. (2.4.69) 

With a, £ 1 , this is always satisfied. For M 2 > 2/(ya,) we require the coefficient of 8 P, to be greater 
than the sum of the coefficients of &P r - 2 and 8 P l+l , 

2/A/ + y M - yA/0,/2 > yMa x - 1/A/ (2.4.70) 

or 

a x < 2/(yA/ 2 ) + 2/3 . (2.4:71) 


2.0 ANALYSIS SECTION 


54 



STABILITY LIMITS. A1=A2=0 



0.00 0.20 0.40 0.20 0.(0 1.00 1.20 1.40 1.00 1.00 2.00 2.20 2.10 2.00 2.00 1.00 

MRCH 


Fig. 2.4.4 Acceptable values of a^ as a function of Mach number 
based on Eqs. 2.4.60, 2.4.62, 2.4.63, and 2.4.65 
( for K = 1.4 ) . 
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Thus we can have a, = 1 up to M 2 = 6/y or up to M = 2.07 for y = 1.4. Fig. 2.4.5 shows the valid 
range of a, based on these criteria. 

We now consider combinations of <%, a, and Oj. In particular if 

+ «3j + 02 = 1. (2.4.72) 

the interpolation scheme is second order accurate. ( See Appendix A. ) 


For Mach numbers less than 2, a, — 1 is stable. Therefore, for M £ 2, we will choose 


02 = 0 

% + a x = 1 . 

From similar stability analyses to those already given 

oi, < 2/(yA/ 2 ) - 1/3 


should be stable for M £ 2. 


(2.4.73) 


(2.4.74) 


For M > 2, we will choose 


4o = 0 
a \ + ** 1 * 

The stability analysis suggests that acceptable values of a, are 

o, < 0.4 + 3.6/(y M 2 ). 

The stability criteria, Eqs. 2.4.74 and 2.4.76, are shown on Fig. 2.4.6. 


(2.4.75) 


(2.4.76) 
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STABILITY LIMITS. A0=A2=0 



o.oo o.ao o.oo o.to o.to i.oo l.ao i.oo i.io i.oo i.oo a.ao a.«o a. so a. so o.oo 

MACH 


Fig. 2.4.5 Acceptable values of a^ as a function of Mach number 
based on Eqs. 2.4.67, 2.4.68, 2.4.69, and 2.4.71 
( for 'i =1.4). 
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0.10 0. tO 0. BO l.oo 1.30 1.10 1.00 1.00 3.00 


STABILITY LIMITS, RO*fll = 1 . A2 = 0; A0=O.fil*A2=l 



O 


0.00 0.30 0.10 0.00 0*10 1.00 1.30 1.10 1.10 1.00 3.00 3.30 3.10 3.00 3.00 0.00 

MACH 


Fig. 2.4.6 


Stability limits for when = 0 and a Q + aj^ - 1, 
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A set of equations for a u and Oj, which satisfy Eq. 2.4.72 and so give second order accurate in- 
terpolation, and which satisfy Eqs. 2.4.73-2.4.76 so that they satisfy the stability criteria, have been 
selected. These are: 


uo = (0.8/3)(4/A/ 2 - 1) 

for M £ 2 a, = 1 - Uq (2.4.77) 

a 2 = 0 ; 

% = 0 

for M > 2 a, = 4/A i 2 (2.4.78) 

oj = 1 - a,. 

These Mach number dependent formulations for <% t o, and a, are shown in Fig. 2.4.7. These 
equations are tested in Section 3.4 where they are referred to as the M&M formula. 


2.5 TIME STEPS 


A unique feature of this method is the use of different time steps for the continuity and momentum 
equations. Previous workers who have used explicit time marching methods have used the CFL 
condition as a basis for determining allowable time steps which maintain stability. The general 
form of the CFL condition for three dimensional inviscid flow is (21), 


( 5 ') cfl * (-^T + + "TET + 


1 


1 


1 


- 1 


[hxf {by) 2 (5z) 2 


( 2 . 5 . 1 ) 


The same time step is used for both the continuity and momentum equations and typically the 
governing equations are updated simultaneously rather than sequentially as is done with the current 
method. The CFL condition is justified by requiring that the analytical domain of influence lie 
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Fig. 2.4.7 


M & M Mach number dependent values ( Eqs. 2.4.77 
and 2.4.78 ) for the coefficients in Eq . 2.4.57. 
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within the numerical domain of influence (21). The CFL condition can become very restrictive 
when the grid becomes highly streched ( the 8 y dimension is much less than the Sx dimension). 
A highly streched grid is needed in calculations with turbulent boundary layers. 

The Denton code (13) uses a variation of the CFL condition (Eq. 2.5.1). The time steps are based 
upon the following, 

5/ £ FT x - t=== - (2.5.2) 

V^o. 

where FT is a time factor typically between 0.2 and 0.5, 8/ min is the minimum characteristic length 
for a given control volume, and y/y RT Ql is a reference speed of sound based upon the inlet stag- 
nation temperature. The same time step is used in the continuity and momentum equations. The 
time steps can ; however, be different for each control volume if only the steady-state solution is 
desired. If a time accurate solution is desired then the same time step must be used everywhere in 
the flow and it must be based upon the control volume with the most restrictive time step. 

The current method not only uses different time steps for each control volume but it also uses dif- 
ferent time steps for the continuity equation and for the momentum equations. The advantage of 
using different time steps is that for flows with thin boundary layers, the allowable momentum time 
step can be significantly larger than that allowed by the CFL condition. These larger time steps 
allow the boundary layer profiles to change more rapidly and enhances the convergence rate sig- 
nificantly compared with a method which uses the CFL condition. To be able to have a calculation 
method which uses different time steps for different equations, the usual conservative form of the 
momentum equations cannot be used. The form of the momentum equations used for the current 
method is 


p-^- + F • pu m — u( F • pw) = — F • P8 U + F*pFu + F*|iF/ . (2.5.3) 

ot J 


This is in contrast to the straight conservative form of the momentum equations which is 
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— ~ + V'p uu- — V ’ P8n + F’pFy + F' Fu T . 

at J 


(2.5.4) 


By using Eq. 2.5.3 instead of Eq. 2.5.4, the a velocity vector is updated directly from the momen- 
tum error. Eq. 2.5.4 requires that you update (pa) first and then calculate the new velocity from 

u n+ ' = (pu) n+ 'lp n+ ' • (2.5.5) 


The time derivative of (p a) , 


<Kpu) 

dt 


can be expanded, using the chain rule, as 


d(P«) _ j. ..3p 

I — P I — • 14 I 


dt 


dt 


dt 


(2.5.6) 


The term dpldt is also calculated from the continuity error. So to correctly update (pa) over one 
time step using the conservative form, Eq. 2.5.4, the same time step must be used for both the 
continuity and momentum equations. However, Eq. 2.5.3 does not have that restriction. It should 
be noted that as the solution approaches a steady-state, (u( V • pa)) goes to zero, therefore Eq. 2.5.3 
then reduces to Eq. 2.5.4. In the current method, the expressions that are used to determine the 
allowable time steps are , for the momentum equations 





(2.5.7) 


and for continuity, 




2 /? 7 ]- htfn 


5/„ 


(5*) 2 (8y) 2 


+ I 


RT 5jc 


I + I 


v eff 

RTdy 


(2.5.8) 


where 5 t m is the momentum time step, 8t c is the continuity time step and v tff is an effective y- 
component of velocity. The effective velocity, v 9ff , is determined from the following equation, 

_ I mass flux through north face I + I mass flux through south face | c 

y eff (25 - 9) 
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where the north and south faces for a control volume are shown in Fig. 2.5.1. These time steps 
are typically reduced by a factor of 2 to reflect the uncertainty of the non-linear nature of the 
equations and the non-uniformity of the grid. During the initial transients, a factor of up to 4 may 
be needed. These reductions in the time steps will be given the symbol TIMEF. In the following 
two subsections, the logic behind these time steps will be presented. The time steps can be derived 
based upon (1) the differential, or (2) the discretized form of the governing equations. 


DIFFERENTIAL FORM 


The x-momentum equation can be written in non-conservative form as 


M + p up- + pv^ = - 

dt dx dy 3x ^ 


(2.5.10) 


assuming that viscous stresses can be represented by the simplified form shown. For stability we 
require that the coefficient of the unsteady velocity be greater than the sum of the absolute value 
of the coefficients which involve derivatives of the u-velocity. Let us rewrite Eq. 2.5.10 in terms 
of changes 5 u, 8t m , and the dimensions 5x, and, 5y, 


8 u 
8t~ 


, 8 u 


5 u = _ 5P 


8u 


i U U . UM _ Ut . UU 

+ + " -ST + 


(2.5.11) 


Now let us apply the above stability requirement to Eq. 2.5.1 1, and the result is 


Tfc > |£| ♦ 4*1 + '^ 


pv 


A L 


(2.5.12) 


or in terms of 8t„, 






+ 


2u 

P(5y) 2 


(2.5.13) 
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The v-velocity used in Eq. 2.5.13 should be the effective velocity which was introduced in Eq. 2.5.9. 


The time step used for the continuity equation is determined using a similar procedure, but because 
of the sequential nature of the current method, extra care must be taken to quantify the interaction 
between the continuity and momentum equations. In the current method the pressure is updated 
because of the continuity error and this new pressure is used in the momentum equation. This 
pressure change must not be too large so as to upset the stability of the method. The continuity 
equation can be written in differential form as 

it + i£!i + iei » o . (2.5.14) 

dt dx dy v ' 


Expanding the derivatives using the chain rule, Eq. 2.5.14 becomes 


dp + du + u dp + dv + v dp _ ^ 

dt P dx dx " dy dy 


(2.5.15) 


What we want to do now is to make each derivative in Eq. 2.5.15 in terms of a common variable. 
Because we are using a pressure update method we will choose pressure as our common variable. 
For our analysis we will assume that the fluid is an ideal gas and that temperature changes are small 
in comparison to pressure changes over one iteration. That means that a change in density, Sp, 
from one iteration to the next can be related to the change in pressure using 

5p = • (2.5.16) 


Rewriting Eq. 2.5.15 in terms of finite changes 8p,8«,8v, and 8/ e and finite dimensions 
8x and 8j> , 


5p 

5'c 




Substituting Eq. 2.5.16 into Eq. 2.5.17, we get 


(2.5.17) 
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The momentum equation must be used to relate a given change in u,8m , to a change in pressure, 
5 P. If the momentum equation, Eq. 2.5.11, is reduced to its steady form it becomes, 


_..8u , 8 m _ __ 5 P j_ 8 m 

+ P'V - -5F + • 


(Syr 


(2.5.19) 


Solving for 8m, and using absolute values to give conservative results, 


8m = 


5P 

P 


l'w' + ^ + 

J p(8y) 


(2.5.20) 


then using Eq. 2.5.13, Eq. 2.5.20 becomes, 


8 u - 


/ 8P \ g» 

_ (— )K, 


OX 


(2.5.21) 


Similarly we can obtain 


/ 8 P 


(2.5.22) 


Now we can substitute Eq. 2.5.21 and Eq. 2.5.22 into Eq. 2.5.18 and the resulting expression is 

1 8/> P . 8^8'm , . _m_ 5/> , P , 5p8t m , , v 5P _ n ™ 

/?T 8* c ”8x l _ p6x 1 /?7 “5F lyP ” pfcy* £7 "5)T * ( • • 3) 


Using the same stability requirement as was used for the momentum equations, we get 


1 ' 1 l-St-l + l-TT^I + I- 8 '" 


RT 5 t c 


(8 xY 


RT Sx 


( 5 / 


V efT 


RTdy 


(2.5.24) 
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which when solved for 8 t c becomes 


8 t c £ 


tf7]| 


6 /. 


(Sx ) 2 


-I + I 


5 /„ 


(8>-) 2 


■I + I 


RTbx 


I + I 


RThy 


11 


(2.5.25) 


DISCRETIZED FORM 


Another way of determining the time step limitations for the continuity and momentum equations 
is to analyse the discretized form of the steady momentum and continuity equations. The mo- 
mentum equations will be analysed first to determine the momentum time steps. 

MOMENTUM TIME STEPS: Integrating the steady flow x-momentum equation 

V • p nu = — + V ' p Vu (2.5.26) 

cx 

over the control volume shown in Fig. 2.5.2 and using Gauss' theorem to obtain area integrals we 
obtain 


^jupu'dd = JJ/M* + fyVu'dd • (2.5.27) 

Discretizing this equation for the control volume shown in Fig. 2.5.2 and neglecting derivatives w'ith 
respect to x in the viscous terms, we obtain 

)ean,U + ('"“k’est.V + ("“Onorrt.W + (™)sou,h,V + (^/+U “ 

- W . . 0 (2 .5 , 8) 

The momentum fluxes, (mu) wttt and {mu) tgst can be represented as, 
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(mu) west = m^ at x u u 


(2.5.29) 


and 

= "W x U 1+IJ • (2.5.30) 

Assuming that the grid spacing is uniform in the y-direction, the momentum fluxes ,(riiu) MIlrt and 
{mu) mrdt are respectively, 

(™ u )south = 0-25 x m sout , | x (u ItJ + u l+ u + 1 ) (2.5.31) 

and 

("*“) north = 0-25 x m north x (u;^ + «/+i^ + «/,/+! + u /+i,/+i). (2.5.32) 

Combining Eqs. 2.5.28, 2.5.29, 2.5.30, 2.5.31, and 2.5.32, we get that 
"W x U I+\J + "W x «// + °’ 25 X "WfA x (“/,/ + U /+V + u U+\ + M /+U+ l) 


+ 0.25 X m south X (u /y + u I+lJ + + u /+1/ _,) + (F/ +1>y - 

- = tftfs . (2.5.33) 

For the converged solution, the right hand side ( RHS ) of Eq. 2.5.33 will be zero. For intermediate 
solutions, the right hand side of Eq. 2.5.33 is equal to the momentum error. 

The stability requirement to be applied to the discretized momentum equation ( Eq. 2.5.33 ) is that 
the change in velocity for node ( 1+ 14 ), 5 u,+ tJ , for one time step be less than the momentum 
error divided by the sum of the coefficients associated with the centerpoint node ( I + 1J ) , in other 
words, 
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* ( momentum error ) 

OW/_j_ j ^ _ ■ 

^ ( sum of the coefficients of node I + 1 ,J ) 


(2.5.34) 


The coefficients associated with the velocity u /+ XJ from Eq. 2.5.33 are 

E( Coefficients ) = m east + 0.25m soulh + 0.25 m„ orth + 2p-|*- . (2.5.35) 

To reflect the uncertainty of the signs of the fluxes through the transverse faces, the transverse 
contributions to the coefficients will be combined into a single transverse mass flux ,m r , defined 
as 

m T — 0.5[| m south | + I m norlh I ] . (2.5.36) 

With this, Eq. 2.5.35 becomes 

Z( Coefficients ) = th east + Q.5m T + 2p-|~- . (2.5.37) 

The mass flux through the east face, m„„, in terms of the fluid properties at the node point (I + 1 J) 
is, • 

™east = P/+ \J U l+W • (2.5.38) 

Therefore Eq. 2.5.37 can be written as, 

L( Coefficients ) = p I+u u l+u 8y + 0.50 m T + 2p-^- . (2.5.39) 


To put the stability requirements of Eq. 2.5.34 in terms of time marching terminology, note that 
from the above stability requirement that 


momentum error 


( sum of the coefficients of node I + IJ ) 
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and from time marching terminology that 


5u, 


7 + 1 ,/ 


_ ( momentum error ) St„ 
P KbZ 


(2.5.41) 


U 


Combining Eqs. 2.5.40 and 2.5.41, we get for a time marching method that 


5/ m £ 


P Vol,j 


( sum of the coefficients of node I + 1 J ) 


(2.5.42) 


Substituting the value for the sum of the coefficients from Eq. 2.5.39 into Eq. 2.5.42, we get that 


5/ m £ 


pS-xSy 


puby + 0.5m T + 


(2.5.43) 


where Vol tJ = 5jc8j; has been used. Simplifying Eq. 2.5.43, we get 


8 'm * 


u + 0-S0^r . 2u 

■5F "pW ?hy i 


(2.5.44) 


Then substituting in the previously defined v tff> ( Eq. 2.5.9 ), Eq. 2.5.44 becomes 


8 £ 


0.5v. 


~5x 




stL 


p8y 2 


(2.5.45) 


which can be compared with Eq. 2.5.13 derived from the differential form of the governing 
equations and is 





(2.5.46) 
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CONTINUITY TIME STEP: The continuity time step can also be calculated based upon the 
discretized form of the momentum and continuity equations. The basic stability requirement used 
in obtaining the continuity time step will be that over one time step that the change in mass flow 
rate not be greater than the error in mass flow rate as defined in section 2.4. Because the degree 
of complexity would be large to analyze the full form of the discretized continuity and momentum 
equations for a generalized control volume, a simplified analysis will be used here which analyzes 
the x and y- direction contributions to the continuity time step separately. The principle of 
superposition will then be used to develop a complete expression for the continuity time step. 

In the PRESSURE INTERPOLATION SECTION ( 2.4 ), the stability characteristics of various 
pressure interpolation schemes were discussed. The stability of these various pressure interpolation 
schemes were analyzed using the discretized form of the continuity equation and the x-momentum 
equation in a 1-D example ( see Fig. 2.5.3). One of the stability requirements introduced in this 
analysis was that the error in mass flow rate and the change in mass flow rate as defined in that 
section have the same sign. 

In the pressure interpolation stability analysis, the discretized form of the steady continuity equation 
was represented by Eq. 2.4.11, which is 

P|+1^/+1® M /+1 “ p /^/5w/ + w /+l / */+l5p/+ 1 ~ = ™error t I (2.5.47) 

where the changes 5( ) are differences between the correct solution and the current solution. The 
left hand side of this equation may also be used to evaluate the changes in net mass flow rate in 
one time step if the changes 5( ) are changes in properties in one time step. Over one time step, 
the changes in pressure are determined from the continuity equation and the changes in velocity 
are determined from the momentum equation. Changes in density are found through the ideal gas 
equation of state. For the one-dimensional example shown in Fig. 2.5.3, if we are using the pres- 
sure update time marching method, along with an upwind effective density, and if we assume that 
only control volume ( 1J ) has a mass flow rate error, we can say that 
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8P, = m errorJ bt c RT)Vol I 


(2.5.48) 


8P ,+ , - 0 (2.5.49) 

5p/ = P t - i/RT = 0 (2.5.50) 

and 

8p, + i = 8Pj/RT = m error j&t c /Vol[ . (2.5.51) 

The changes in velocity are determined from the discretized form of the momentum equation. The 
resulting expressions for the change in velocity, 8u, were found to be 

8u /+ , = (8P, - 8 P I+ i)p /+ ^ X/ = [rtierrotjS'c P/+ ”b Xl (2.5.52) 

and 

Su, - <«/>,- , - «*,) t, e (2.5.53) 

having assumed that at the beginning of a time step that the momentum equation is balanced except 
for the continuity error. Now if Eqs. 2.5.50, 2.5.51, 2.5.52, and 2.5.53 are substituted into Eq. 
2.5.47, we get 


^ a r R.T i 

P/+ \ A I+ 1 I m error J ° l c y Q ^ J p /+ , $ X{ 


RT 


8t„ 


PjAj j m e r TO , j 5t c 1 p/8jr ; _ 


8r c 

+ Uj+ ]/!/+ terror, I m change 


(2.5.54) 
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For our one-dimensional example, 


Vol, S A fix, (2.5.55) 

A l+t ^A, . (2.5.56) 

5x/_ , £ 5x, (2.5.57) 

and 

P/+1 = P / • (2.5.58) 

Using Eqs. 2.5.55, 2.5.56, 2.5.57, and 2.5.58 and rearranging, Eq. 2.5.54 becomes, 

1 • (2.5.59) 

The upwinded density method ( q, = 0, a, = 0, a, = 0 ) has been used in our analysis because 
it results in the most conservative time step. Eq. 2.5.59 is used in the PRESSURE INTERPO- 
LATION SECTION to see under what conditions the error in mass flow rate has the same sign 
as the change in mass flow rate. Now we want to And an equation which tests under what condi- 
tions the change in mass flow rate is less than the error in mass flow rate, 

™ change ,/ < ™ error J • (2.5.60) 

Substituting Eq. 2.5.60, into Eq. 2.5.59 we get 

1.0 * • (2-5.61) 

Solving Eq. 2.5.61 for the continuity time step, 5 t c , we get 
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1 


(2.5.62) 


bt c <: 


RT\ + u 


[ 


5jc 


2 RTbx 


] 


When the differential form of the governing equations was used, we got 


5 t e ^ 


Rn ±HL + + |_iL_| + I-^Z-I] 

8x 2 by 2 RT5x RTb y 


(2.5.63) 


where the terms in Eq. 2.5.62 can be recognized as the x-momentum contributions to Eq. 2.5.63. 


THE CFL CONDITION AND THE DENSITY UPDATE METHOD: When the density update 
time marching calculation procedure is used by Denton ( 13 ) f the same time step, the CFL con- 
dition, is used in both the continuity and momentum equations. It will be shown using the stability 
criteria for the continuity time step , < m „ roT ), that if the same time step is used in both the 

momentum and continuity equations, Eq. 2.5.62 gives a result close to the CFL condition. It 
should be noted though that the density update method is not necessarily restricted to using the 
same time step. 

Using an analysis similar to that used in the previous section, it can be shown that if the density 
update method is used along with the downwind effective pressure, the resulting one dimensional 
relationship between changes in mass flow rate and errors in mass flow rate is 

l 2 + § x RT ]™ error *I ^ c ^ change (2.5.64) 

which is identical to Eq. 2.5.59. Now we want to determine what the limiting stable time step is 
if we let 5 t m = 8t c and < m error . The resulting expression for the time step is 

5 ' * ~ U + ^ + 8/?71 • (2-5.65) 
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For our one-dimensional example, the CFL condition would set the time step as 


8 trFL* ~ - • (2.5.66) 

CFL {u + c) 

Fig. 2.5.4 shows a comparison of the (time steps/Sx) calculated using Eqs. 2.5.65 and 2.5.66 for a 
Mach number range from 0.01 to 2.0. Over this range, the time steps calculated using these two 
equations are within ± 20 %. At Mach numbers greater than 0.4, the CFL condition gives a more 
conservative time step than Eq. 2.5.65. Fig. 2.5.5 shows a comparison of the time steps calculated 
using Eqs. 2.5.65 and 2.5.66 with the different time steps for the continuity and momentum 
equations calculated using Eqs. 2.5.7 and 2.5.8 over the same Mach number range. 

TRANSVERSE CONTRIBUTIONS TO THE CONTINUITY TIME STEP: When mass im- 
balances are caused by incorrect mass fluxes through the transverse faces, the pressure changes from 
the continuity error act through the y-momentum equation to correct these mass imbalances. An 
analysis, similar to that just shown for the 1-D example, will be used to show how transverse fluxes 
affect the continuity time step. 

Fig. 2.5.6 shows four control volumes of identical dimensions at a fixed axial location. All errors 
in transverse fluxes are assumed to be zero except between contol volumes ( 14 ) and ( 14-1 ). To 
the right of the control volumes, in Fig. 2.5.6, is a plot of the static pressures at the upstream nodes 
for the control volumes, before and after one time step. The static pressures at all nodes are as- 
sumed to be initially the same ( represented as triangles in Fig. 2.5.6 ). 

The transverse fluxes shown in Fig. 2.5.6 will cause the pressure (represented as squares in Fig. 
2.5.6) to increase at node (14-1) and to decrease at node (14) through the mass flow rate errors for 
control volumes (14*1) and (14) respectively. These changes in pressure at nodes (14) and (14-1) 
will induce pressure gradients in the transverse direction (represented by the dashed line in Fig. 
2.5.6) which will cause the mass imbalance to correct itself. The pressure changes ; however, must 
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TIME STEPS 




Fig. 2.5.5 Time Steps Comparison 
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A INITIAL PRESSURE AT NODES 
□ FINAL PRESSURE AT NODES 
O FINAL INTERPOLATED PRESSURES 


Fig. 2.5.6 Simplified Geometry For Determination of 
Y-Momentum Equation Contribution to Time 
Step 
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not be so great so as to violate our stability requirement that the change in mass flow rate be less 
than the error in mass flow rate. 


Several simplifying assumptions are made throughout this part of the analysis to reduce its com- 
plexity. Just as in the analysis in the PRESSURE INTERPOLATION SECTION, we will assume 
here that the momentum equation is balanced initially for each control volume except for the 
continuity error. In Section 2.4, the x-momentum equation reduced to Eq. 2.4.24 

5 «,+ , = (5 P, - hP I+ ,)5f/(p /+ ,8*,) (2.5.67) 

when the momentum equation is balanced initially except for the continuity error. For the case 
where the transverse fluxes cause the imbalance, the equivalent y-momentum equation is 

8v/+i,/ = 0-25 * (8 P ij + 8P/ + u + §Pij- 1 + SPj+i^-i) 

-0.25 x (5 PlJ + 5 P l+lJ + 8P IJ+X + 5 P /+u+1 ) ■ . (2.5.68) 

1P/+ ijtyij) 

Eq. 2.5.68 can be simplified for this uniform grid to 

5v /+u = [0.25 x (5 P IJ . l + 5 - 0.25 x (5 P lJ+ , + 8P l+u+ 0 \ pJ ^jE yJ J < 2 - 5 - 69 ) 


For the case that we are investigating, 

8 Pij - , = - 5 t c RT/Vol, (2.5.70) 

8P /+U _, = 0 (2.5.71) 

5P V+I = 0 (2.5.72) 

SJW+i-0. (2-5.73) 
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Therefore, in this case, Eq. 2.5.69 becomes 


8v /+u ( m„ TOr 0.25 Vo ^ ) p/+i " 6j;/ 


(2.5.74) 


The same procedure can be applied to control volumes ( IJ-1 ) and ( 14 + 1 ) and the resulting 
expressions for changes in v are 


Rv =< — m x 0 25 x 8t w 

5v /+U-l ( m error 0- y^ ) p /+ ^ 


(2.5.75) 


and 


8 y /+ 1^+ 1 - i™error,IJ x 0.25 x 


ML 5^ 

Ko/ v J P/+V+1^/ ' 


(2.5.76) 


The change in mass flow rate due to the changes in v for control volume ( 14 ) is determined from 
the continuity equation using the y-contribution equivalent to Eq. 2.4. 1 1 in Section 2.4, which 
would be 


change (p/45v) n0 rth + iy^P)north (Vj48p) J(?u ^ * (2.5.77) 


The change in mass flow rate for one iteration for our uniform grid is therefore, 


m change 


p x 0.25 x 5x r x (5v /f/ 4 5v /+u 4 5 v /y+1 4 8v /+If</+1 ) 

— P X 0.25 X §X[ X (5V[j 4- 8v /+Jj y 4* Sv^.j 4 5V; +It y_j) 

4 V X 0.25 X X (5p/^ 4 Spj+J^y 4 Sp;^ +1 + 5P/+l,/+l) 

— V X 0.25 X &Xj X (Spjj 4 5p/ +lf y 4 Sp/^y^j 4 8p /+lj y_j) 


(2.5.78) 


assuming mean values for the density and v-velocity. The density changes for one iteration are 


5p /f y ™ error, [j§t c / Volj 


(2.5.79) 


Sp/y- 1 “ “ ™error,l t fit c IV°l l • 


(2.5.80) 


2.0 ANALYSIS SECTION 


82 



Substituting Eqs. 2.5.74, 2.5.75 2J.76, 2.5.79, and 2.5.80 into Eq. 2.5.78, we get 


change IP terror y^j ) P( "terror foe j/„/ ) 1 oSvT - X 


Vol, " P Wl 


&‘c 


v bx/ ( 0-25m error j i j ) • 


(2.5.81) 


If Eq. 2.5.81 is rearranged, we get 


_ °- 25 x "terror RT , 0.5ht m hx, ^ V 5x, , 

m change - y Q , i ^ dt 1 ' 


*7 


(2.5.82) 


For the control volumes used in this analysis, 


Volt = 8 jc /5 y t 


(2.5.83) 


Applying the stability criteria, m ehmt , < m,„„, and using Eq. 2.5.83, Eq. 2.5.82 becomes, 


1.0 * 8 r e * 71 -^^ + 4^-1 ■ 


Syf 


RTZx, 


(2.5.84) 


When Eq. 2.5.84 is solved for the continuity time step, 8t e , we get 




(2.5.85) 


The mass flow rate error through the transverse face of control volume ( IJ ) causes the pressure 
to change at node ( I,J ). This pressure change will also act through the x-momentum equation in 
exactly the same manner as previously outlined in Eqs. 2.5.47-2.5.63. When both the x and y 
contributions to the continuity time step ( Eqs. 2.5.62 and 2.5.85 ) are combined together we get 
that 
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1 


a 


(2.5.86) 


5 t e £ 


RT[^f 

5x 2 


1255f m 

5/ 


+ 


u 

RT8x 


+ 


Q.25v 

RT&y 


To insure that the time step we choose will result in a stable calculation procedure, the time step 
that has been calculated using the differential form of the governing equations ( Eq. 2.5.25 ) and the 
time step derived from the descietized form of the governing equations ( Eq. 2.5.86 ) have been 
combined together and symmetry has been invoked and the resulting equation is Eq. 2.5.8 pre- 
sented earlier, 


5 t c z 


2 RT[^~ 

8x 2 


K 

5 / 




(2.5.87) 


RT8x RTZy 


2,6 BOUNDARY CONDITIONS AND INITIAL 
GUESS 


BOUNDARY CONDITIONS 


Along the upstream boundary, the total temperature, total pressure, and v-velocity are specified for 
inviscid flow. Along the downstream boundary the static pressure is specified. For viscous flow, 
at the upstream boundary, the total temperature, freestream total pressure, inlet boundary layer 
velocity profile, and flow angle are specified. 
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Pressures along the solid boundaries are determined from linear extrapolation. There is no mass 
flux across element faces which coincide with the solid boundaries. For viscous flow, the values 
of the x-component and y-component of velocity are set equal to zero at solid walls. 

For flow through cascades, the additional boundary condition of periodicity must be considered. 
Fig. 2.6.1 shows a two dimensional projection of a typical grid system up to the leading edge of a 
cascade blade. Note that a grid point is not located along the periodic boundary when the new 
control volumes are used. The computational domain extends from the lower periodic boundary 
to the upper periodic boundary. The missing calculation points outside the computational domain 
are replaced by the corresponding points adjacent to the other periodic boundary. 


INITIAL GUESS 


The initial guess for the inviscid part of the flow field is determined from a 1-D inviscid solution. 
A boundary layer is then added along the wall using a constant ratio of boundary layer thickness 
to duct height throughout the duct. The velocity profile used in the boundary layer is the inlet 
velocity profile. For certain geometries, an estimate of the blockage effect of the boundary layer is 
used to specify an effective geometry for the calculation of the initial solution. 


2.7 CALCULA TION OF VISCOUS FORCES 


The momentum equation for unsteady flow in differential form is 


d£U_ + 
dt 


Vp ‘uu= — VP+ V'\iVu + 


pVu T + 


V 




(2.7.1) 
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The normal stresses associated with second coefficient of viscosity, X, will be neglected in this 
analysis and in any subsequent calculations. We will be concerned in this section with how the 
viscous terms (F*pF& and V'pVy 7 ) are evaluated using the non-orthogonal, physical mesh 
system that is incorporated in the present method. 

In a control volume analysis of a flow field, we are interested in the actual forces which act upon 
the control surfaces and the components of these force in the coordinate directions ( x,y, and z ) 
rather than the derivatives of the shear stresses as seen in Eq. 2.7.1. To transform the governing 
equations from differential form (Eq. 2.7. 1) into an integral form, we use Gauss' theorem. In terms 
of some arbitrary vector, g>. Gauss' theorem says, 

JJJ V ’ SB dVol = JJa> • dd . (2.7.2) 

The viscous terms ( V • p Vy and V • p VyT) in Eq. 2.7.1, are converted from differential form to 
control volume form using Eq. 2.7.2 and the result is, 

JjJ( V • p Vu + V • yivifidVol = JJ(p Vu-dd + \i~Fu T ‘ dd) . (2:7.3) 

The current two-dimensional scheme uses control volumes which are made up of four straight line 
segments ( see Fig. 2.7.1 ). The surface integral in Eq. 2.7.3 is simplified into a summation over 
the four sides of the control volume. The integral over the surface in Eq. 2.7.3 can therefore be 
represented as 

JJ(p Vu'dd +\LVu T 'dd) = 

I Uk * H + Ak ‘ H Pa 7 ) • (2-7.4) 

K= 1 

Each area vector in Eq. 2.7.4 is assigned a magnitude equal to the area of the face in question and 
a direction which points in the outward normal direction. The evaluation of Eq. 2.7.4 results in the 
net viscous forces in the x and y directions for a control volume. 
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Fig. 2.7.1 Typical Control Volume 
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For laminar flow, the absolute viscosity, p, is used in Eq. 2.7.4 to evaluate the shear forces. For 
turbulent flow, an effective viscosity is used in Eq. 2.7.4. The effective viscosity that is used is the 
sum of the absolute viscosity and an apparent eddy viscosity. The current method uses a Prandtl 
mixing length model to evaluate this eddy viscosity. This mixing length model is outlined in Table 
2.7.1. 

The details of this mixing length calculation will be presented later in this section. 

The velocity gradients , V u L , needed in Eq. 2.7.4 can be determined within a non-orthogonal grid 
by using, 


y u - x ^ U L + &K x d u L + Hi x Dj 3u l 

L £>r(Dj * Ok) dl Q r (Qj X el k ) dJ £,•(£, * R K ) dK 

where 22/, 22/, and 22* are directional vectors along the grid directions (IJ, and K ), see Appendix 
B 0 For the two-dimensional case, duJdK = 0 and 22* is a vector spanning the height of the duct 
( in the direction 22/ x 22/). Two typical two dimensional control volumes are shown in Fig. 2.7.2. 
The directional vectors ( 22/ and 22/ ) are identified in Fig. 2.7.2. The magnitudes of the vectors are 
dependent upon the grid spacing as can be seen in Fig. 2.7.2. The directional vectors that are shown 
in Fig. 2.7.2 would be used to calculate the velocity gradients applicable to the boundary common 
to both of the control volumes. 

The derivatives of the velocities in Eq. 2.7.5 are taken with respect to the grid indices, in other 
words, 


” ( u l)/+i ~ ( u l)i (2.7.6) 

where (u L ) ; is the velocity at the beginning of vector 22/ and (w L ) /+ j is the velocity at the end of 
vector 22/. For the geometries and boundary conditions investigated in the present work, the gra- 
dients in properties in the I direction are much smaller than gradients in the J direction. As a 
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Table 2.7.1 Prandtl Mixing Length Model 


Vta ■ + m 


(i) 


= 


# u 

dy 


( 2 ) 


L is the smaller of 


0.08 times the width of the boundary layer 

/ 

0.41 times the distance to the nearest wall 


Van Driest Correction 


L = 0.41T( 1 “ exp[ ~ YyfinWvto < 3 ) 


Near Wall Correction 


^eJT = yj^i + 53 


14 ) 
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consequence, only the J derivative contributions to Eq. 2.7.5 will be considered in the present work. 
This is equivalent to using the thin-layer Navier-Stokes formulation. 

FORCES ON THE SOUTH FACE OF A CONTROL VOLUME 

In the actual calculations, for a non-uniform , non-orthogonal grid, the directional vectors are 
slightly different than those shown in Fig. 2.7.2. For the south face of a control volume, the area 
and directional vectors actually used are shown in Fig. 2.7.3. The directional vector J2j, and the 
velocity change are evaluated using the downstream nodal values because their use strengthens the 
centerpoint coefficient of the matrix of unknown variables. The directional vector J2 h is located 
midway between the four nodes rather than at the boundary surface. This is because the viscosity 
and velocity gradient used in calculating the shear forces on the south face are evaluated midway 
between the two nodes. The shear stress is known to vary less through the boundary layer than the 
velocity gradient or the mixing length squared. Therefore, it is preferable to calculate the shear 
stress using a velocity gradient and mixing length midway between the grid points in the J-direction 
and then assign the resulting shear force to the face of the control volume between the points. The 
upper wall and lower wall control volumes are shown in Figs. 2.7.4 and 2.7.5 , respectively with the 
directional vectors identified. The shear stresses are evaluated midway between the wall and the 
near wall point and then the shear forces are assigned to the wall surface. The effective viscosity 
used to evaluate these wall shear forces is a combination of the laminar and turbulent viscosities 
given by 


V-eff = VVM + Ht) • (2-7.7) 

This relationship is used only at the wall and has been shown ( Ref. 22 ) to allow a good calculation 
of wall shear stress with a near wall point further away from the wall than is typically required. 

Symmetry is used to calculate the forces on the north face of a control volume, in other words, 

FnorthJrJ ~~ Fs 0U thjj + 1 (2.7.8) 


2.0 ANALYSIS SECTION 


92 




Fig. 2.7.3 Two Typical Control Volumes With Directional 
Vectors Identified For South Face 
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where F north and F indl are the shear forces on the north and south faces of the control volume re- 
spectively. 

FORCES ON THE WEST FACES OF A CONTROL VOLUME 


The forces on the west face of the control volume are calculated by scalar multiplying an interpo- 
lated shear stress with the west face area vector ( see Eq. 2.7.4 ). Refering to Fig. 2.7.6, for a node 
point (14) which is located on the west face of control volume (14). first the velocity gradient and 
viscosity are calculated at the western edge of surfaces A and B which are midway between the node 
points ( 14 ) and ( 14 * 1 ) and node points ( 14 + 1 ) and ( 14 ). respectively. The west sides of surfaces 
A and B are also the locations where the shear stresses are calculated for the north and south faces 
of the control volume (1-14)- The product of the velocity gradient and viscosity at node ( 14 ) is 
determined by linearly interpolating using the following interpolation formula, 


<PV " <Pa + 


<fc/2) 


y 2 2 ' 


(*Ps ^/l) 


(2.7.9) 


where q> A and <p B are the products of the velocity gradient and the viscosity at the west sides of 
surfaces A and B respectively. Once these interpolated values have been calculated, Eq. 2.7.4 can 
be used to calculate the components of the shear forces on the west face. The shear stress is in- 
terpolated to node 14 rather than interpolating the velocity gradient since the shear stress varies less 
through the boundary layer than the velocity gradient. Similarly for the east face we have 

^eastJJ ~ ~ F\vestJ+ \J • ( 2 . 7 . 10 ) 


CALCULATION OF MIXING LENGTH AND VISCOSITY 

When the mixing length model of Table 2.7.1 is used in the 0.4Iy region, a distance normal to the 
wall, "y”, must be determined to calculate the mixing length. For a flat wall, the distance to the 
wall from a point is, of course, measured along a line perpendicular to the wall. However, with the 
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Fig. 2.7.7 Distances Normal to the Wall 
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current grid system, grid lines and lines orthogonal to the wall are not coincident ( see Fig. 2.7.7 ). 
Fig. 2.7.7 shows 5 adjacent control volumes. The normal distance (L) between the north and south 
faces of a control volume is equal to the volume of the control volume divided by the distance (S) 
between the east and west faces of the control volume measured in the I grid direction. The total 
normal distance from the wall to the node point on the west face of a control volume ( ) is 

equal to the sum of all previous normal distances between that control volume and the nearest wall 
plus one half of the normal distance for that control volume. We can represent that distance as 

^west = + 0.5 x Ln (2.7.11) 

k=\ 

where n is the number of control volumes from the wall. 

The total normal distance needed for calculating the mixing length for the south face of a control 
volume ( L south ) is equal to the sum of all previous normal distances to the previous node ( 
L wnt ^ x ) plus the equivalent normal distance half way between the two grid points in question. 
The procedure used to calculate this length is shown in Fig. 2.7.8. We can represent this total 
distance as 


Lsouth j i - 1 + 


(A,-l +Ln) 


(2.7.12) 


To determine the mixing length in the outer part of the boundary layer, the boundary layer thick- 
ness measured normal to the wall must also be determined. The lengths that were calculated above 
can also be used in calculating the boundary layer thickness. The edge of the boundary layer is 
determined by using the magnitude of the normalized local total pressure gradient as a measure of 
its location. The normalized total pressure gradient used here is 


normalized total pressure gradient 


A P, l_ 

'A r 0.5p e u\ 


(2.7.13) 
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where 'A y m is the normal distance between two grid points. The local freestream values of density 
and velocity are used in Eq. 2.7.13. The grid line that is considered the freestream is input into the 
method. Starting from the freestream the normalized total pressure gradient ( Eq. 2.7.13 ) is cal- 
culated between grid points and is compared with a characteristic normalized total pressure gradient. 
The edge of the boundary layer is defined as where the magnitude of the calculated gradient is equal 
to the characteristic gradient. The exact distance is determined using, interpolation and is based 
upon the distances normal to the wall. The characteristic normalized total pressure gradient is de- 
termined by 


^t,char 1 _ 1 

~W~ 0.5p A CL 


(2.7.14) 


where CL is some characteristic length of the flow field, typically the duct width. 


Once the mixing length has been calculated, the apparent eddy viscosity is calculated using Eq. 1 
in Table 2.7.1. The velocity gradient used in the eddy viscosity calculation is determined using the 
following expression, 


"dif _ /dL* Fuj i + , 4* P~v .2 

dy r ui ur. 


(2.7.15) 


where A is the area vector on the south side of the control volume. This formula gives the mag- 
nitude of the total velocity gradient reflecting the non-orthogonality of the grid. 


2.0 ANALYSIS SECTION 


101 



2.8 MULTI-VOLUME METHOD FOR PRESSURE 


CHANGES 


Preliminary calculations of boundary layer flows using the density update method (with non- 
uniform grid spacing) resulted in solutions which became unstable after only a small number of it- 
eration steps. After a detailed investigation of the nature of this instability, the cause could be 
attributed to effects resulting from the large aspect ratio of the control volumes near the wall and 
the large variation of properties in the boundary layer. The first contributing factor was the use of 
different time steps for each control volume and for each equation. Near the walls where the con- 
trol volumes were very thin and the velocities were low, the continuity time steps that were used 
for calculating the changes in density were small. The ratio of the momentum and continuity time 
steps there was also very large. The pressure in the outer part of the boundary layer was changing 
much more rapidly than it was near the wall. This induced large transverse pressure gradients and 
then large transverse velocities followed resulting in an unstable calculation procedure. It was felt 
that to stabilize this calculation procedure, the changes in pressure through the boundary layer must 
be coupled in some manner and that the changes in pressure be only dependent on the continuity 
error and not on both the density change through the continuity error and the temperature change 
through the momentum error and its resulting velocity change. We wanted to minimize transverse 
pressure gradients in the intermediate solution to enhance stability. 

The above realizations resulted in two changes. One change altered the way that the continuity 
error is used to update the flow properties (see Section 2.3 DISTRIBUTION OF PROPERTIES). 
Previously (12), errors in continuity were used to update the density at the node points. The 
pressure was then calculated from the equation of state. An alternative procedure has been devel- 
oped which updates the pressure directly from the continuity error. The density is then evaluated 
using the equation of state. 
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The second change is to group control volumes in the boundary layer to form a larger global con- 
trol volume. The continuity error is calculated for this global control volume and changes in 
pressure are assigned equally to each of the upstream nodes for each control volume making up the 
global control volume. Initially the global control volume extends from the wall to the edge of the 
boundary layer. Then the global control volume is made successively smaller towards the wall. 
This is shown schematically in Fig. 2.8.1. The entire pressure change for one iteration at each node 
within the multi-volume region is determined by adding together all the pressure changes assigned 
to that node . If a non-uniform grid is used, the continuity time step for each global control volume 
may based upon the continuity time step for the largest control volume within the global control 
volume. A non-uniform grid is preferred for this multi-volume method so that pressure changes 
for successively smaller global control volumes become smaller. 

The multi-volume method propagates pressure changes rapidly through the boundary layer and 
minimizes large tranverse pressure gradients in the intermediate solution. The above changes allow 
the calculation of boundary layer flows where the control volumes near the wall can have aspect 
ratios (length/height) of over 1000. Without the multi- volume method these aspect ratios could 
not be used. The stability and convergence properties of various multi-volume approaches are in- 
vestigated in section 3.5. 


2.9 TRANVERSE UPWIND DIFFERENCING 


When the control volumes become long and thin near the wall of the duct, the fluxes through the 
top and bottom faces of the control volume become more significant in comparison to the fluxes 
through the streamwise faces. Because the nodes of the control volumes are located in the middle 
of these streamwise faces, the predominant flow direction must be in the streamwise direction for 
the discretization method used here to properly reflect the convective nature of the flow. When the 


2.0 ANALYSIS SECTION 


103 






fluxes in the transverse direction become significant, the solution procedure may become unstable 
because the diagonal terms in the coefficient matrix become smaller as the transverse fluxes increase. 
This is because the velocities at the comers of the nodes are determined from interpolation. To 
strengthen the diagonal dominance of the coefficient matrix, the momentum fluxes through the 
transverse faces may be calculated using interpolated velocities upstream in the transverse direction 
rather than the actual interpolated values. These velocities are multiplied by the mass fluxes 
through the sides of the control volumes to get the total momentum flux. The direction and 
magnitude of the upwinding is determined from a criterion based upon the magnitude of the flux 
ratio, FRATIO. FRATIO is the ratio of the mass flux through the transverse face divided by the 
mass flux through the streamwise face. A criterion for upwinding is derived next. 

The criterion for upwinding is based upon an analysis of the convective terms in the momentum 
equations and how the discretization affects the dominance of the centerpoint coefficient. The 
momentum equations for unsteady flow are 

~=~ + V'p yy = — V ' P8n + V'pVy + V'pVu 1 ^ . (2.9.1) 

of J 

When Gauss' theorem is applied to the steady convective terms of Eq. 2.9. 1 to transform the gov- 
erning equations from differential to integral form, we get 

JJJ V • pyyd(Vot) = jjpyy • 44 . (2.9.2) 

The convective terms in Eq. 2.9.2 can be rewritten as 

B = p4=- + V • py y — y V • pa . (2.9.3) 

ot 

When Gauss' theorem is applied to the steady portion of Eq. 2.9.3, we get 

C = JJJ( V • py u — y V • pg) dVol = * 44 ~ uJjpM • 44 (2.9.4) 
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where U is an average velocity for the control volume. The second term on the right hand side of 
Eq. 2.9.4 can be recognized as the continuity error contribution to the momentum equation. We 
use Eq. 2.9.4 in the current method, to evaluate the momentum fluxes. 

If our computational domain is descretized into finite quadrilateral control volumes, the convective 
momentum fluxes can be identified as those associated with the north, south, east, and west faces 
of the control volume ( see Fig. 2.9.1 ). These momentum fluxes can be further classified into those 
associated with the x-momentum equation ({rini) mra ,,(mu), ouA ,(mu ) w „ tt and ( mu) ml ) and those asso- 
ciated with the y-momentum equation ((Am')„, rt ,(mv) reBrt ,(mv) w „ t , and(mv) M „) . For example, the 
momentum flux contribution to the x-momentum equation through the north face of the control 
volume is 


(rr* u)north - pa • 44 (2.9.5) 

where A is an outward normal vector. A typical control volume with the x-momentum fluxes 
identified is shown in Fig. 2.9.1. 

The mass fluxes through the north, south, east, and west faces of a control volume are 
n%, ort K m, oytK m,„ t% and ^ respectively. The mass flux through the north face of the control 
volume out of the control volume is 


"Wrt “ (pU'dA (2.9.6) 

where positive mass flux is defined here as a flux directed outward from the control volume. The 
continuity error, m 9rron for a control volume (IJ) is 

^error ” — ™northJJ ~ ™ southJJ ^westJJ ~ ™ eastJrf * (2.9.7) 

From symmetry, we can note that 
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m north JJ fn southJ r J+\ 


(2.9.8) 


and 

^ eastJrJ ~ ™westj+ \J • (2.9.9) 

Substituting Eqs. 2.9.8 and 2.9.9 into Eq. 2.9.7, we get 

™error JJ ~ ^houthJJ+A ~~ ^houthJJ — ™westJrf • (2.9.10) 

From the momentum balance, the convective contributions to the x-momentum equation ( see 
Eq. 2.9.4 ) for a finite control volume are 


C (fhu )northJff south JJ ^^eastJJ + {^E^erro^IJ (2.9.11) 

where (w £ m tffror ) w is the continuity error contribution to the momentum error. The corrections to 
the momentum error, {u E m error ) u and {v E m trror ) u , due to the continuity error are 

(terror) U = - «JJp U‘ dd (2.9.12) 


and 


tfpu-dA (2.9.13) 

where u E and v £ are effective velocities for a control volume and U and v are average velocities for 
a control volume. The effective velocities may be different from the average velocities to improve 
the stability of the method. 

The momentum flux (mu)„ r o, > s > 
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(PU * d) north J,J u northJJ 


{™H) n orthJJ ^btorthJjUnorthJJ 

“ — nhouthJJ+xUnorthJJ (2.9.14) 

Using a four point interpolation for 

(rtw) north Jfl “ — **houthJJ + 1 x / 

+ W/+i/) x (1 “ gpj) + ( u /,y+i + u /+i,/+i) x 4 ?at1 (2.9.15) 

where 0 ^ g N £ 1. Similarly, 

= ™southJrJ x + w /+l,/-l) x U “£$) + ( w /,/ + U !+\J) x £sl (2.9.16) 

i™ u )westjj = "W,/,/ x «/,y (2.9.17) 

{jhlfyeastJJ *“ “ (^ w )vves7,/4* \J ^Wr/,7+1,/ X U I+ \J (2.9.18) 

where g N and g s are interpolation parameters for the north and south faces of the control volume 
respectively. The interpolation parameters , g H and g s , are usually calculated from the position of 
the control volume face relative to the two node points adjacent to that face (see 2.2 CONTROL 
VOLUMES). We will call this the geometric interpolation parameter, g^. However when the 
transverse fluxes become large enough in comparison to the streamwise fluxes, the centerpoint co- 
efficient will not be dominant if these geometric interpolation parameters are used. Other more 
stable interpolation parameters must be determined. 

Now substituting Eqs. 2.9.10, 2.9.15, 2.9.16, 2.9.17, and 2.9.18 into Eq. 2.9.11, and grouping coef- 
ficients with common velocities, we get 
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C - U l+\J * l ™southJ,J + 1 x 0.5 x (1 - g N ) + rhswthJJ x 0-5 x gs ~ r KiesiJ+ 1/1 
+ U u * [ “ "houlhJJ+l X 0.5 X (1 - gff ) + x 0.5 x fo + "Vert,/,/! 

+ U IJ+ 1 x I ~ ril south,lJ + 1 x 0-5 x g N ] 

+ u I+\J+\ x ( - ri houth,IJ+ 1 x 0-5 x 
+ u U-\ x [™south,V x 0.5 x (1 - fo )J 
+ u l+\J-\ x l " 1 * outhJJ x 0-5 x (1 - g 5 )] 

+ Mg x [ + r ^;outh,IJ+ 1 — ™southJJ ™westj+ \ J ~ ™west,IJ I • (2.9.19) 

The value of u £ in Eq. 2.9.19 will depend upon the sign of the continuity error and its influence on 
the strength of the centerpoint coefficient. 

The last term in Eq. 2.9.19 is 

+ Mg x | + m S0Uth j fJ+ j — m south Ir/ + fewest j + l J ~ f^westjA 

= + Mg x [m error \ . (2.9.20) 

If the continuity error, m, rror , is negative we will let 

Mg = u u (2.9.21) 

which reduces the magnitude of the non-centerpoint coefficients. If the continuity error, m „ nt , is 
positive we will let 

“£=“/+!/ (2.9.22) 
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which will strengthen the centerpoint coefficient. Both these decisions enhance the stability of the 
method. It can be shown for forward flow and with this choice of the effective velocity, u E , that 
the coefficient of u u is always negative if either = 0 or m,outhjj+ \ = 0- As will be seen, 

the contribution of node IJ to the sum of the non-centerpoint coefficients will not need to be 
considered further in this analysis. Also because of this choice of u £ , the magnitude of the mass flux, 
does not detrimentally affect the stability of the scheme in forward flow. 

To arrive at interpolation parameters which will probably result in a stable calculation procedure, 
we will investigate the stability of two simple examples. These examples assume 1) that the mass 
flux through the north face of the control volume is zero and 2) that the mass flux through the 
south face of the control volume is zero. Finally, more conservative interpolation parameters will 
be determined which will take into account the possibility of arbitrary combinations of north and 
south fluxes. This will be accomplished by using only 1/2 of the flux, in the calculation 

of the interpolation parameters which satisfy the stability criteria. 

CASE # 1 

For this case, we have assumed that the mass flux through the north face of the control volume is 
zero . This test case is shown schematically in Fig. 2.9.2. When the mass balance is applied we 
get, 


™errorJJ ™southJJ ^\vestJJ ™westJ+\rJ • (2.9.23) 

If the mass flux through the south face is positive and the continuity error is negative because of 
this flux, then Eq. 2.9.19 becomes 

C * «/+U X l°' 5 * Ss * "'southJJ ~ MwestJ+W + « error 1 
+ Ujj X [0.5 * gs x "houthJJ + "Wr/,/,/1 
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+ U;,/-! x |0.5 X ihsouthjj x (1 - g 5 )] 


+ u l+\J-\ * 10-5 x outhjj x 0 ~ &s)l • 

All non-centerpoint coefficients are negative; therefore, the requirement for stability is 

0.5 x g s x rn zouth l j ~ ril w est,l+\ r J + terror * 0 
0.5 x g s x m, 0UthJrl - ~ rh southJJ > 0 . 

This will be satisfied for any value of g r 

For positive rn, 0 uth.u and a corresponding negative continuity error, Eq. 2.9.19 becomes 
C = U l+\J x x 8s x ™south,lJ ~ ™west,I+\,jl 
■*" u ij x [0.5 x gs X ™south,IJ rawest ,IJ ™error\ 

+ U[j — ! x [0.5 x m south jj x (1 - g s )] 

+ M /+l,/-i x (0-5 x nhouthjj x (1 ~ £s)l • 

The requirement for stability is, therefore, 

0.5 x gs * "houthjj ~ ™west,i+\j ^ ™southjj x (1 “ Ss) • 

In terms of FRATIO, fK 0 uu,jjl f K n t,i+\j< and gs, Eq. 2.9.28 becomes 

> 1 + FRATIO 
8s 1.5 FRATIO ' 


CASE # 2 


(2.9.24) 

(2.9.25) 

(2.9.26) 


(2.9.27) 

(2.9.28) 

(2.9.29) 
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For this case, we have assumed that the mass flux through the south face of the control volume is 
zero. When the mass balance is applied we get, 

^errorJJ ~ ^soulhJJ + 1 — ™west,IJ ™west,I + 1 J • (2.9.30) 

If the mass flux through the north face is positive, in other words, if m t0UtKlJ + l < 0 and is 

negative then Eq. 2.9.19 becomes 

C ~ u l+\r! x I “ ™southJJ+\ x 0-5 x (l — gfi/) - ™westj+ \j\ 


x f ** houthJJ+l x x (1 Sn) revest, IJ ^errorl 


+ U u+i x l g N x 0.5 x m south jj+ x J 

+ u l+\J+\ x [ - ^ x 0.5 x m S0Uth lfJ+x J . (2.9.31) 

The requirement for stability is, 

— r *houth,I r J + 1 x 0.5 x (1 - g N ) — fn wes tj+ X j ^ g/y x m s 0 uth,Ij+\ • (2.9.32) 

In terms of FRATIO and g N , we get 


< 1 + 0,5 x FRATIO 
gN 1.5 FRATIO 


(2.9.33) 


.vhere in this case, 


FRATIO = ntsou,h ' , ’ J+ 1 . 

f * t westJ+ \J 


(2.9.34) 


If the mass flux through the north face is negative, in other words, if m I0Uth%u+i > 0 and rn, mtj j is 
positive, then Eq. 2.9.19 becomes 
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^ U l+\J x 1 ^southj J + 1 x 0.5 x (1 <?/v) ™westj + 1 J ^errorl 


+ U U x l~ m southJJ+ , X 0.5 * (1 - g N ) + m westj j\ 


+ U U+\ x I &N x 0-5 X ^southJrl+x] 

+ u l+ i/+l x [ — gff x 0.5 x m south j r j + j) . (2.9.35) 

The requirement for stability is, 

“ ™south,ij + 1 x 0.5 x (l - gN ) - m westJ+w + m error > 0 (2.9.36) 

— ™southjj+\ x 0-5 x (1 — gpj) — rhyygffjj + ^houthjj + 1 ^ 0 . (2.9.37) 

This result is always satisfied for all values of g N . 

Since the north face of one control volume is the south face of another control volume, 
3 ^,+, = g NJ . The most conservative value for the interpolation parameter is chosen from Eqs. 
2.9.29 and 2.9.33 so that we use a consistent interpolation scheme for a given face. If the geometric 
interpolation parameter , g INT , is more conservative than those specified in Eqs. 2.9.29 and 2.9.33, 
then it should , of course, be used. 

ANY COMBINATION OF NORTH AND SOUTH FLUXES 


If only 1/2 the mass flux, is used in the determination of the interpolation parameters, 

g N and g s , it can easily be seen ( looking at Eqs. 2.9.28 and 2.9.32 ) that the corresponding in- 
terpolation parameters become more conservative and are 


^ 0.5 ± FRATIO 
8s 1.5 FRATIO 


(2.9.38) 
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£ 0.5 ± 0.5 FRATIO 
8n l.SFRATIO 


(2.9.39) 


This is equivalent to using twice the FRATIO in Eqs. 2.9.29 and 2.9.33. The criteria for the in- 
terpolation parameters, g s an&g N , are summarized in Table 2.9.1. The above interpolation pa- 
rameters are not assured to provide stable discretization for all possible combinations of mass fluxes 
; however, these interpolation parameters have resulted in stable solutions for all the test cases in- 
vestigated in this dissertation. However, some of the solutions were unstable before this upwinding 
was used. For all the test cases in this dissertation, upwinding was only needed in the transient 
solution and at the steady state the geometric interpolation parameters were used. 


2.10 ENERGY EQUA TION 


For most of the calculations which are to be used as test cases in the present work, the assumption 
of constant total temperature will be a sufficient representation of the energy equation in the flow 
fields. By assuming constant total temperature, the computations are less expensive to run and the 
computer storage requirements are less. An assumption of constant total temperature is satisfactory 
because of these reasons: 

1. an adiabatic wall is assumed in the calculations. 

2. no work is done on the fluid at the solid boundaries. 

3. the Mach numbers in the flow fields that are investigated are low enough that total temperature 
gradients within the boundary layer are small. 
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Table 2.9.1 Upwinded Interpolation Parameters 


SOUTH FACE 


If 

m , . . < 0 

south, i,j 

no restriction on g g 

If 

m . . . > 0 

south,i,j ' 

0.5 + FRATIO 


1 . 5FRATI0 


where FRATIO = m south,i,j 

* 

m . , , 

west , i+l , j 


NORTH FACE 


If m u • • . i < 0 
south, i, j+1 


If m . . . t . >0 

south, l, j+1 


g -10.5 + 0,5 FRATIO 
n 

1.5 FRATIO 
no restriction on g 

n 


m 


where FRATIO = south , i , j+1 

m . , , 

west , l+l , j 


- noting that for consistency g for control volume j 

n 

must equal for control volume j+1 

reflecting symmetry and the geometric interpolation 

parameters, the conservative value is chosen 
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4. the turbulent Prandtl number of air is approximately 0.9. For a Prandlt number of 1.0, the 
assumption of constant total temperature is valid, therefore a solution with a Prandtl number 
of 0.9 should not deviate greatly from the constant total temperature assumption. 

However, it was felt that the methodology required to include the energy equation in the governing 
equations should be developed and simple test cases should be used to demonstrate its success. 
Later in this dissertation, four test cases will be presented which use the full energy equation in their 
calculations. These test cases will be 

1. turbulent boundary layer flow in an adverse pressure gradient with an inlet freestream Mach 
number of 0.55. 

2. flat plate turbulent boundary layer with a freestream Mach number of 0.95. 

3. Sajben's diffuser (23) including the full energy equation. 

4. flat plate turbulent boundary layer with a freestream Mach number of 2.8. 

Unless otherwise noted, the laminar Prandtl number for the calculations will be .73 and the tur- 
bulent Prandtl number will be 0.9. For the moment, the integral form of the energy equation in- 
corporated into this finite volume method will be presented. 
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(2.10.3) 


— + V -E# = - — + F • pe^ 


5(pe«) ,, , de, „ 

a ~~ + P • (pe,)a = p-^- + pu • Ve, 


(2.10.4) 


and then expanding the right hand side of Eq. 2.10.4, we get, 


de, _ de, _ , _ 

— + p u' Ve,= p— + V • pue, - <?,( F • pu) 


(2.10.5) 


The procedure just outlined is identical to what was done to the unsteady and convective terms in 
the momentum equations (see section 2. 1). 


The heat flux vector can be represented as 


5 - — k VT 


( 2 . 10 . 6 ) 


Substituting Eqs. 2.10.3-2.10.6 into Eq. 2.10.1 , we get 


de, _ 

— - - V • pi/e, + 


F • pue t + e,( V • pu) — F • ( — k VT) 


+ F • [u ■ (p Vu + p Vm t ) J - V 'Pu • 


(2.10.7) 


Applying Gauss' theorem to convert to the integral form of the energy equation we get 


x 8 Vol = - JJpue, • c/4 + <?,JJpu • dd. ~ JJ — k VT • dd 


+ JJ[u • (p Fu + p Fu 7 )) • dd ~ JJFu • dd 
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where e t is an average value for the control volume. As with the momentum equations, Eq. 2.10.8 
has a term, y • dA , which removes the continuity error contribution to the energy error. 

To enhance stability, the energy equation of Eq. 2.10.8 is not used from the initial solution. The 
first 500 iterations of the calculation use the assumption of constant total temperature. After 500 
iterations, the energy equation is added. Initially there are large errors in continuity and momentum 
and these large errors act through the energy equation to cause errors in the total energy for a 
control volume. This interaction is destabilizing therefore we let the calculation proceed until 
continuity and momentum are reasonably satisfied before the energy equation is added. 

The appropriate boundary conditions to be used with the energy equation in our calculations are: 

1. Adiabatic wall, meaning that the heat flux through the solid walls is zero. 

2. There is no work done on the control volumes at the solid boundary due to viscous forces. 
This is because the velocity at the wall is zero. 

Because the energy equation has the same general form as the momentum equations, approximately 
the same time step is used in the energy equation as is used in the momentum equations. More 
will be said about this later in this section. The addition of the energy equation does destabilize the 
solution somewhat when compared with the assumption of constant total temperature. For this 
reason, the time step used for the energy equation is usually reduced by a factor of 2 compared to 
the momentum equation time step. Transverse upwind differencing is also used for the specific 
total energy . e v in the energy equation in exactly the same manner as it was used for velocities in 
the momentum equations. 

An alternative form of the energy equation has also been used and will be derived here. This al- 
ternative form has enhanced convergence properties when compared with the above formulation. 
Briefly, the energy equation is reformulated so that changes in total enthalpy ,/t,, are calculated over 
one time step rather than changes in total energy ,e, , which was done previously (Eq. 2.10.7). Also, 
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the energy error due to the continuity error and the momentum error are removed using this alter- 
native formulation. 


The total enthalpy can be defined in terms of the total energy and the static temperature by starting 
with the definition of total enthalpy 


h t = 




(2.10.9) 


and for an ideal gas, 

t~ rt 


( 2 . 10 . 10 ) 


so that 


h t =e, + RT . (2.10.11) 

Taking the derivative of both sides of Eq. 2.10.1 1 with respect to time , t , and multiplying by the 
density, we get 




+ p R 


dT 

dt 


( 2 . 10 . 12 ) 


The static temperature T can be represented in terms of the total enthalpy and the absolute velocity 
as 



therefore 


8T _ \ dh t V dV 

dt C p dt C p dt 


Therefore Eq. 2.10.12 can be rewritten as 


(2.10.13) 


(2.10.14) 
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Now if we solve for the time derivative of the total enthalpy, we get 



= pyifL _ P* v dV 


dt 


dt 


(2.10.16) 


where y is the ratio of specific heats and V is the magnitude of the velocity vector. The second term 

on the right hand side of Eq. 2.10. 16 is the momentum error contribution to the total enthalpy er- 
dc 

vox. Substituting p-r*- from Eq. 2.10.7 into Eq. 2.10.16, we get 
at 



— V ' pye t + e,( F • pa) — F • ( — k VT) 


+ F • [a • (p Vu + p Fa 7 )] - 


.CL dt 


(2.10.17) 


This equation may be rewritten, by using ft, = e, + P/p, as 



— V * pu/t t + h t ( V • pa) — V • ( — k VT) 


+ F-[a-(pFa + pFa r )i - -£-F- P a - • 


(2.10.18) 


By noting that V • pu — dp/dt, the last two terms of Eq. 2.10.18 are of the form, 


M = 




(2.10.19) 


At the steady state, Eq. 2.10.18 becomes 

0 = y [ - V • pafy - V • ( — k VT) + V • [a • (n Fa + p Fa 7 )]] . (2. 10.20) 
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Therefore, we may arbitrarily alter the variables 1 and m and the steady form of the energy equation, 
Eq. 2.10.20, will be obtained for converged solutions. The transient behavior of A, is improved in 
the calculation procedure by choosing 1 = m = 0 so that the transient energy equation used for 
each time step is, 



V • p yh t + h t { V • pu) — V • ( — k VT) 


+ V •[«' (p Vu + p Fu r )l . 


( 2 . 10 . 21 ) 


By setting 1 = m = 0, we remove the influence of the continuity and momentum errors to the 
energy equation during transients in the solution. Applying Gauss' theorem to Eq. 2.10.21 we get 
the integral form of the energy equation, 



8 Vo/ = y 


— fjpyh t ’ d£ + A,JJpy * dd. ~ JJ “ k VT • dd 


+ JJ[a • (p V U + |i V U T ) 1 * dA (2. 10.22) 

where h t and P are average values for the control volume. The time step used in Eq. 2.10.22 must 
be reduced by a factor of y compared with that used in Eq. 2.10.8. This reduced time step is needed 
to ensure the stability of the new formulation since the coefficient of the steady terms in the energy 
equation have increased by a factor of y in the new formulation (see section 2.5). 


It was found to be important in the implementation of the energy equation, that the heat flux term 
and the work done due to shear forces term were evaluated consistently. This means that the work 
due to shear forces should be evaluated midway between the node points and not at the control 
volume faces because the heat flux calculation by the nature of its discretization is valid at the 
midpoint between the node points. 


An interesting observation can be made from this form of the energy equation about boundary layer 
flow over a flate plate with a Prandlt number of 1. For steady flow, Eq. 2.10.22 becomes 
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0 = — JJp uh, • dA ~ JJ — kVT ' dA 


+ JJ[« * (p Vu + P Fy 7 )] ' dA 


(2.10.23) 


In addition, the viscous terms are simplified using the boundary layer approximation by noting that 


U‘ (p Vu + p Vu T ) s up-|j- = P-^- 


(2.10.24) 


and therefore 


dip'll 


JJ[« • (p Vu + p Fa 7 )] • dA = JJp-^^4 


(2.10.25) 


For an ideal gas with constant specific heats, 


/l 


(2.10.26) 


and the Prandtl number is defined as, 




(2.10.27) 


therefore, 


- - kVT-dA = JJ-^- Vh-dA . (2. 10.28) 

Also if gradients in enthalpy in the x-direction are small compared to gradients in the y-direction, 
then 

(2.10.29) 

Substituting Eqs. 2.10.24, 2.10.25, 2.10.28, and 2.10.29 into Eq. 2.10.23, we get 
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(2.10.30) 


0 = - jjpuh, • dd + JJ-£- 2LdA y + JJp-^ • 

The last two terms of Eq. 2.10.30 can be combined if the Prandtl number is taken as 1. The re- 
sulting equation is 


0= - JJp Uh t -dd + . 

If the total enthalpy is constant in the flow, Vh, = 0, therefore 

Sl^y " 0 


(2.10.31) 


(2.10.32) 


and 


jjpu/t, • dd = ~ h t j$pu‘dd = 0 


(2.10.33) 


at the steady state. Therefore the specification of constant total enthalpy in the flow field is a sol- 
ution to the energy equation for flow over a flate plate when the Prandtl number is 1. It should 
be noted also that the term, is a diffusion term in terms of the total enthalpy ,/t,. The 

energy equation and the momentum equation are therefore similar enough to justify the use of 
approximately the same time steps for both the energy and momentum equations. 
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3.0 RESULTS 


3.1 FLA T PLA TE BOUNDARY LA YER 


A laminar boundary layer was calculated in a constant height duct. The boundary layer thickness 
at the inlet was 15 % of the duct height. The freestream Mach number was 0.43. The inlet velocity 
profile was the Blasius profile. The absolute viscosity was 0.0 1 kg/m s. The duct height as 44 mm 
and the duct length was 112 mm. The geometry and the grid are shown in Fig. 3.1.1. The 
Reynolds number based upon x varies from 5070 to 11840 along the duct. The duct is 17 inlet 
boundary layer thicknesses long. Fig. 3.1.2 shows that the development of the velocity profile 
compares very well with that predicted by theory. 
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3.2 INVISCID CALCULA TIONS OF SAJBEN'S 
DIFFUSER 


Sajben's diffuser (23) will be used here as a test case to illustrate the difference in the results obtained 
using the effective pressure and effective density methods (see section 2.4). Remember that an ef- 
fective pressure or an effective density is needed to stabilize the solution procedure for the density 
update and the pressure update methods, respectively. The calculations in this section were made 
using a three dimensional finite- volume code derived from the Denton finite volume code (13) now 
in use at NASA Lewis Research Center. This new code uses the current contiol volumes (see 
section 2.2) in a three dimensional configuration. 

The geometry and grid used in the calculations are shown in Fig. 3.2.1. There were 34 axial grid 
points and 10 equally spaced radial grid points. The current calculations are made essentially two 
dimensional by inputing the coordinates of the diffuser at a very large radius ( 900 m ) in x-r co- 
ordinates. The calculations begin at x/h = -3.6 and end at x/h = 7.9, where h is the throat height. 
The inlet total pressure is 135 kPa and the inlet total temperature is 300 K. The exit static pressure 
is 108 kPa. This gives a P 0Xi JP ttM9t = 0.800. With these conditions, one dimensional isentropic flow 
gives a shock with an upstream Mach number of 1.495 at the location marked in Fig. 3.2.1. 
Multigridding (13) is used to improve the convergence speed. 

Fig. 3.2.2 shows a comparison of the effective pressure and the thermodynamic pressure for this test 
case when the the density update method was used. In regions where the flows varies smoothly, 
the effective pressure and thermodynamic pressure are essentially equal but in regions of the flow 
where there are large gradients in properties they can be substantially different ( up to 10 % ). Fig. 
3.2.3 shows a comparison of the total pressure through the diffuser calculated using the 
thermodynamic pressure and the effective pressure. The total pressure calculated from the effective 
pressure is in much better agreement with the theoretical solution; however, current practice is to 
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print out the total pressure calculated from the thermodynamic pressure. The exit P U jJP t jniwt IS 
calculated to be .931 for both cases and it agrees well with the one dimensional analytical value of 
.9304. 

Fig. 3.2.4 compares the bottom flat wall static pressures for calculations of transonic flow in 
Sajben's diffuser using the effective pressure method and the effective density method with the same 
boundary conditions as specified above. The pressure shown in Fig. 3.2.4 from the density update 
method is the effective pressure. Fig. 3.2.5 compares the total pressure for these two cases. The 
effective density method gives a much more uniform total pressure upstream and downstream of 
the shock; there are no overshoots in total pressure when the effective density method is used 

In summary, much better total pressure distributions through shocks are obtained when the in- 
terpolated effective pressure, needed to stabilize the density update solution procedure, is used to 
calculate the total pressure. This simple change largely eliminates the undershoot in total pressure 
downstream of a shock. Overshoots and undershoots in total pressure can then be further reduced 
by a factor of 10 by adopting the effective density method rather than the effective pressure method. 


3.3 THE INFLUENCE OF TRANSVERSE 
SMOOTHING ON A STEP PROFILE IN A 
STRAIGHT DUCT 


Transverse smoothing is required in Denton's method with the control volumes shown in Fig. 
2.1.2 because there are more grid points across the duct (unknowns) than there are control volumes 
(equations). Smoothing formulae are used to add non-physical 'extra equations". Two forms of 
transverse smoothing are used in the Denton code; these are linear smoothing, described in Table 
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Comparison of Total Pressures Calculated Using The Effective 
Pressure and Effective Density Methods for Sajben's Diffuser 



3.3.1, and non-linear smoothing, described in Table 3.3.2 and Fig. 3.3.1. This transverse smoothing 
of properties cause numerical viscosity to be introduced into the solution when large gradients in 
the properties are seen across the duct. 

For calculations where the gradients across the duct are small, little numerical viscosity would be 
expected. However, in the boundary layer region of a viscous flow the gradients in properties can 
be large and this smoothing could cause large numerical viscosity. As a severe test case, calculations 
were made for a step profile in inlet properties in a straight duct. Both the control volumes used 
by Denton and the control volumes used in the present work were used in the calculations. The 
geometry can be seen in Fig. 3.3.2. A step inlet profile of total pressure is specified. The total 
pressure at the centerline is 135 kPa and the total pressure is reduced to 120 kPa 
[PuJP,*. nuriin. = 0.889 ) at the sides ( see Fig. 3.3.3 ). The exit static pressure in the duct is 108 
kPa ( 0.8 P,, cmt „ liJU ). 

Fig. 3.3.4 shows Mach number profiles at three axial locations along the duct for the case where 
linear smoothing was used ( with SF = 0.02 ). The inlet step profile ( x = 0.0 m ) is quickly altered 
into a parabolic type profile ( x = 4.0 m ). This parabolic profile then changes relatively little until 
the end of the duct ( x = 21.0 ). Fig. 3.3.5 presents the total pressure distribution along the duct. 
The step profile causes an almost step change in the total pressure at the beginning of the duct and 
then the total pressure decreases as in a viscous flow. Fig. 3.3.6 compares the Mach number profiles 
at the end of the duct for calculations using linear smoothing ( SF = 0.02 ) and non-linear 
smoothing ( SF = 0.02 ). Non-linear smoothing did not improve the profile. 

Additional calculations were made using the same boundary conditions as above but using the new 
control volumes and no smoothing. Fig. 3.3.7 compares the inlet Mach number and exit Mach 
number profiles for this test case. The improvement over the previous results is dramatic. The total 
pressure distribution has also improved especially along the centerline of the duct as can be seen in 
Fig. 3.3.8. These results show conclusively that the numerical scheme used to calculate flows with 
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Table 3.3.1 Linear Smoothing of Flow Properties 


D(J) = (l.-SF) * D(J) + SF*( D(J+1)+D(J-1)) 

2 . 


the variable D at node J is smoothed using this 
equation. The smoothing factor is SF, typically 
0.01 - 0.02. The variables are updated and then 
smoothed. The variables that are smoothed are 
, PV„, {) , and j>e. 
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Table 3.3.2 Procedure for Non-Linear Smoothing 


1) an average value of a property D is determined from 
the neighboring nodes using linear interpolation. 

AVG(J) (see Eq. 1 Fig. 3.3.1 ) 

2) the difference between the actual and average value of 

a property D at a node is determined and assigned the 
variable name CURVE (J). ( see Eq. 2 Fig. 3.3.1 ) 

3) a variable SCURVE is determined from the average of the 
variable CURVE from the neighboring nodes. 

( see Eq. 3 Fig. 3.3.1 ) 

4) the variable D at node I is smoothed using equation 4 
in Fig. 3.3.1. The smoothing factor is SF, 
typically 0.01 - 0.02. 

5) this non-linear smoothing procedure results in no 
smoothing added to linearly or parabolically varying 
properties . 
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D(J+1) 



1) AVG(J)=FD(J)*D(J+1)+FU(J)*D(J-1) 

2) CURVE (J)=D(J)-AVG(J) 

3 ) S CURVE ( J ) = FU ( J ) *CURVE ( J - 1 ) +FD ( J ) *C UR VE ( J+ 1 ) 

4) D(J)= (1 .-SF)*D(J)+SF*(AVG(J)+SCURVE(J)) 

Fig. 3.3.1 Non-Linear Smoothing 
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Fig. 3.3.3 Inlet Total Pressure Profile 
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Fig. 3. 3. A Mach Number Profiles At Different Axial 
Locations With Linear Smoothing 
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Total Pressure Along Grid Lino 
With Linear Smoothing 



MflCH NUMBER 



Fig. 3.3.6 Exit Mach Number Profile Comparison 
Non-Linear and Linear Smoothing 
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No Smoothing - New Control Volumes 


large transverse gradients in properties , like those seen in turbulent boundary layers, must not have 
smoothing of properties in the transverse direction. 


3.4 1-D COMPUTATIONAL TESTS OF SHOCK 
CAPTURING USING PRESSURE INTERPOLATION 
FORMULAE TO CALCULA TE EFFECTIVE 
DENSITY 

DENTON'S 1-D NOZZLE FOR TESTING SHOCK CAPTURING 

Denton ( 12 ) has tested shock capturing with his finite- volume method in a convergent-divergent 
nozzle ( see Fig. 3.4.1 ) designed to produce a linear variation of Mach number with distance for 
1-D isentropic flow. The equation for the Mach number variation with distance is 

* = 10. + 45.(A/ - 1) . (3.4.1) 

Denton considered flow which began at x = 1, where the Mach number was 0.8 and ended at 
x = 46, where the Mach number was 1.8. The throat ( M = 1.0 ) was located at x = 10. He used 
three back pressures with / > „„// > Mn ,.» = 0.85, 0.80, and 0.75, respectively. The theoretical 1-D sol- 
utions for these flows are shown in Fig. 3.4.2. The maximum Mach numbers, just upstream of the 
shock , are 1.268, 1.455, and 1.578, respectively. This is a range of shock Mach numbers typical 
of turbomachinery flows. 

These three pressure ratios are used with Denton's 1-D nozzle to test shock capturing with three 
of the pressure interpolation methods discussed in Section 2.4 ( Pressure Interpolation ). 
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PRESSURE INTERPOLATION SCHEMES 


The effective density method used in the current work uses an interpolated approximation for the 
pressure in the evaluation of the density. A general form of the interpolation formula considered 
in this thesis is 


K M -Pt + ^/-M - Pi) + -yVi+1 - Pi- i) + 1 - ^- 2 ) (3-4.2) 

and this is used to evaluate the density as 

• <3A3) 

This general form is a linear combination of a single-point interpolation, P m — P it a 2 - point in- 
terpolation, — /Vj, and a 3 - point interpolation, P l+1 - /V_ 2 . The single-point interpolation, 
of course, really gives the correct perfect gas equation and involves no approximation. 

The coefficients <^, a u and Oj are here taken to be constants or functions of Mach number. Com- 
binations, including individual terms of pairs of terms, for which the sum of the coefficients 

Oq + + 02 = 1 

are second order accurate, as shown in Appendix A. 

Correct Perfect Gas Equation ( a, = 0 , a, = 0, and Oj = 0 ) 

This scheme has the advantage that it involves no interpolation or approximation for the pressure. 
Experience has shown that it is stable for subsonic flow. But the stability analysis of Section 2.4 
shows that for Mach numbers above about 1.2 this scheme becomes unstable. Thus it could not 
be used for the test cases of Denton's 1-D nozzle. 
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These observations about the use of the correct perfect gas equation are in agreement with Denton's 
findings for his scheme B ( 12 ). In that method changes in density were sent to the upstream 
comers of the element, which is equivalent to our sending pressure changes upstream. The method 
'proved stable, without any correction factors or damping, at low Mach numbers but instability 
was found to develop at Mach numbers around unity and above.' 

2- Point Interpolation ( <% = 0 t a x = 1 = 0 ) 

The stability analysis presented in Section 2.4 shows the 2-point scheme 
numbers up to about 2.0. Use of a 2-point scheme or a 3-point scheme 
Denton in his recent ASME and AGARD Lecture Notes ( 13,14 ). 

3- Point Interpolation ( = 0 ,a t = 0 ,Oj = 1 ) 

3-point schemes have been shown in Section 2.4 to be the most conservative ( in lerms of stability 
) of the schemes considered in this thesis. Perhaps for this reason, such a method is used to stabilize 
the current NASA version of the Denton code. Both 2-point and 3-point schemes provide second 
order accuracy for a continuously changing pressure; they give correct interpolation values for linear 
variations in pressure ( assuming equally spaced grid points). 

M&M Mach Number Dependent Interpolation 

The advantages of the three schemes just considered are : 

1 . the accuracy and stability of the perfect gas equation for subsonic flow; 

2. the stability of the 3-point interpolation at Mach numbers greater than 2.0; 

3. the stability and reduced smearing of properties of the 2-point interpolation at supersonic 
Mach numbers up to 2.0. 


to be stable for Mach 
has been suggested by 
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These advantages have been combined in a single Mach number dependent interpolation scheme 
in Section 2.4. 


In this method 


* ^M 2 


~ 1 ) 


a i = 1 - Oq 


M < 2.0 


a 2 = 0 
% ~ 0 


a, = -4r M > 2.0 

' M 2 

02 = 1 — a x (3.4.4) 

The values of the three coefficients are shown graphically in Fig. 3.4.3. Note, once again, that the 
sum of the coefficients is equal to one for all Mach numbers, so that this scheme is also second 
order accurate. For the calculations presented here M was taken as the larger of M on the upstream 
and downstream side of the control volume. 


COMPUTATIONAL TESTS OF THREE PRESSURE INTERPOLATION SCHEMES 

Of the four schemes just considered, three are stable in the Mach number range 1.0 to 2.0. These 
are the 2-point, 3-point, and M&M interpolation methods. In this section, results of shock cap- 
turing with these three methods are presented and compared for Denton's 1-D nozzle. 

The one-dimensional calculations were performed with 46 axial grid points with an equal axial 
spacing between grid points of 1 unit. The inlet Mach number was 0.8. The ratio of specific heats, 
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7, was 1,4 and the gas constant , R , was 287. J/kg K. Calculations were made with P 9Xit IP t M*t °f 
0.85, 0.80, and 0.75. 

The variations of static pressure , Mach number, and total pressure are plotted for each interpo- 
lation scheme using the same scales as for the theoretical solutions which can be seen in Fig. 3.4.2. 
Fig. 3.4.4 shows the results for the 3-point scheme, Fig. 3.4.5 for the 2-point scheme, and Fig. 3.4.6 
for the M&M method. The results from the 3-point and M&M schemes are shown together with 
the theoretical solution on Fig. 3.4.7 for the pressure ratio of 0.80. 

The calculated values of maximum Mach numbers upstream of the shock and total pressure ratios 
are compared with the values from the theoretical 1-D solutions in Table 3.4.1. 


The total pressure ratios across the shocks are well calculated by all three interpolation formulae 
as shown in Table 2b. This is in spite of the fact that the calculated values for the maximum Mach 
number upstream of the shocks are significantly different from the theoretical values. For example, 
at the lowest back pressure, the theoretical Mach number upstream of the shock is 1,578 while the 
3-point interpolation formula gives 1.501, the 2-point formula gives 1.526, and the M&M formula 
gives 1.533. For this case the calculated values of total pressure ratio are 0.9029, compared with 
the theoretical value of 0.9032. In general the M&M formula gives the closest agreement with the 
upstream Mach number while the 3-point formula gives the worst results. Based on the maximum 
calculated upstream Mach number for these cases , the M&M formula would give shock losses 
from i 6 to 41 percent too small, while the 3-poini formula would give values from 27 to 67 percent 


too small. Interestingly the agreeement for shock losses based on maximum upstream Mach 
number improves ( for all three formulae ) as the Mach number increases. However, these results 
show that the peak calculated Mach number should not be used to predict shock losses and that 
the calculated total pressure loss across the shock is accurate to better than 0.1% and it should be 
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Table 3.4.1. Results of Calculations for Denton f 1-D Nozzle, 
Table 3.4.1a: Maximum Mach number upstream of shock 

Interpolation Formula 


exit 

Theoretical 

3-Point 

2 -Point 

Mach Number 

> 

t inlet 




Dependent 

0.85 

1.268 

1. 175 

1. 198 

1.218 

0.80 

1.455 

1.374 

1.395 

1.408 

0.75 

1.578 

1.501 

1.526 

1.533 


Table 3.4.1b: 

Total pressure 

ratio, P 

t 

. /P 

exit t linlet 


Interpolation Formula 


P 

exit 

Theoretical 

3-Point 

2-Point 

Mach Number 

p 

t inlet 




Dependent 

0.85 

0.9846 

0.9845 

0.9845 

0.9845 

0.80 

0.9433 

0.9431 

0.9432 

0.9432 

0.75 

0.9032 

0.9029 

0.9029 

0.9029 
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The results in Figs. 3.4.4, 3.4.5, and 3.4.6 show that the interpolation formulae all act to smooth 
properties upstream of the shocks. The smoothing is most noticeable in the static pressure and 
Mach number distributions, especially with the 3-point interpolation scheme. The 2-point scheme 
gives less smoothing while the M&M formula gives the sharpest and most accurate distributions. 

Both the 3-point and 2-point interpolation schemes give overshoots in static pressure and under- 
shoots in Mach number downstream of the shocks. Only the M&M interpolation formula shows 
no noticeable overshoots and undershoots and this is because it has a better formulation for sub- 
sonic flow; in fact, from Eq. 3.4.4, it can be seen that the M&M formula reduces to the correct 
perfect gas equation for Mach numbers less than 0.918. 

The M&M formula captures the shocks over about four grid points centered around the theoretical 
shock location. This is seen for the pressure ratio of 0.8 in Fig. 3.4.7. In contrast, Fig. 3.4.7 shows 
the 3-point scheme smearing the shock over about ten grid points with the shock displaced slightly 
downstream due to inadequate resolution of the subsonic flow. Once again the 2-point scheme 
gives results intermediate between those of the M&M and 3-point schemes. 


3.5 LAMINAR BOUNDARY LA YER IN TWO 
CONVERGING DUCTS 


A laminar boundary layer is calculated on the curved wall of two converging ducts. The other wall 
is straight and was treated as inviscid in the calculations. In ail the calculations, the inlet boundary 
layer thickness is 5% of the inlet duct height. The inlet velocity profile is the Blasius profile. The 
inlet height of the duct is 44 mm and the exit height is 31 mm. The length of the duct is 180.4 
mm. The absolute viscosity is .001 kg/m s. The inlet freestream Mach number is 0.10. 
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Two geometries are investigated with the basic dimensions given above. For one geometry, the 
curved wall radius is determined using a sine wave formulation. This geometry has a smooth 
transition from the inlet to the exit sections. The second geometry is essentially the same except 
that the radius is not determined using an analytical function and the second derivative of the wall 
radius is discontinuous. Both uniform and non-uniform grids are used in the boundary layer region. 
Figures 3.5.1 through 3.5.3 show the grid and geometry for the 3 arrangements to be investigated. 
The arrangements are 

1. Smooth geometry using non-uniform grid ( Fig. 3.5.1 ) 

2. Smooth geometry using uniform grid ( Fig. 3.5.2 ) 

3. Non-smooth geometry using uniform grid ( Fig. 3.5.3 ) 

A plot of the 2nd derivative of the wall radius as a function of axial position for the non-smooth 
and smooth geometries is seen in Fig. 3.5.4. The location of the discontinuity is also identified in 
Fig. 3.5.3. 

The specific ideas which are to be illustrated using this test case are, 

1. Comparison of skin friction coefficients calculated using the finite volume method with those 
calculated using Thwaites method. 

2. Comparison of skin friction coefficients calculated using uniform and non-uniform grids. 

3. Comparison of skin friction coefficients calculated using smooth and non-smooth geometries. 

4. Investigation of pressure variations across the duct comparing the uniform and non-uniform 
grids results, and non-smooth geometry results. 

5. Investigation of inlet and exit boundary condition specifications. 
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Fig. 3. 5. A Second Derivative of Radius For Converging Duct 
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6. Investigation of the convergence properties of uniform vs non-uniform grids, coarse vs fine 
grids, and different multi-volume approaches. 

The radial distribution of grid points for the coarse, uniform and non-uniform grids and the fine 
non-uniform grid are tabulated in Table 3.5.1. 

For the smooth geometry in Fig. 3.5.1, the skin friction coefficient along the duct was calculated 
using the finite volume method and using the method of Thwaites (25). A coarse non-uniform grid 
was used for the finite volume calculations. Thwaites method is an integral method used to calcu- 
late the development of laminar boundary layers in incompressible flow. Using the wall static 
pressure distribution from the finite volume calculation as a boundary condition , the skin friction 
coefficient was calculated with Thwaites method. Fig. 3.5.5 compares the skin friction coefficients 
obtained from the finite volume calculations with those obtained using Thwaites method. The 
agreement is good for most of the flow ;however, near the exit of the duct there is a discrepancy. 
This discrepancy may be due to the fact that the boundary layer has grown to over 25% of the duct 
height at the exit. 

The skin friction coefficients calculated with the finite volume method using uniform and non- 
uniform grids are shown in Fig. 3.5.6. The agreement is excellent and it is concluded that grid 
spacing does not affect the prediction of the skin friction coefficient for this laminar flow. 

The skin friction coefficients for the flow through the smooth and non-smooth geometries are 
shown in Fig. 3.5.7. Superficially, the geometries in Figs. 3.5.2 and 3.5.3 look almost identical ; 
however, for the non-smooth geometry the second derivative of the lower wall radius of curvature 
is discontinuous (see Fig. 3.5.4) causing the solution to be erratic at the discontinuity. To ensure 
smooth calculations, geometries input into this program should have continuous second derivatives. 

The cause for the rapid increase in the skin friction coefficient at the discontinuity can be seen in 
Fig. 3.5.8. Fig. 3.5.8 is a plot of a transverse pressure coefficient which is defined as 
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Table 3.5.1 Transverse Distribution of Grid Points 
for Laminar Calculations 



Uniform 

Grid 

Coarse 

Non-Uniform 

Grid 

Coarse 

Non-Uniform 

Grid 

Fine 


y/h 

v/h 

y/h 

PI 

.005 
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.0004 

P2 

.015 

.0071 

.00115 

P3 

.025 

.0142 

.00225 

P4 

.035 

.0283 

.0040 

P5 

.045 

.0563 

.0075 

P6 

.060 

.1125 

.0150 

P7 

.085 

.225 

.0300 

P8 

.125 

.400 

.0575 

P9 

.200 

.625 

.1125 

P10 

.375 

.875 

.225 

Pll 

.625 

- 

.400 

P12 

.875 

- 

.625 
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Fig. 3.5.7 Skin Friction Coefficient Smooth Geometry Vs Non-Smooth Ceometry 
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(3.5.1) 


C Pt 


p_ZIj± 

0.5 p« 2 


where p is the local static pressure at an axial location, p f , is the freest ream static pressure at this 
axial location , and 0.5 pm 2 is the dynamic pressure at the inlet to the nozzle. The different curves 
in Fig. 3.5.8 represent the pressure coefficient at every other grid line from the wall to the freestream. 
Because of curvature, the pressure gradient across the duct first rises and then falls. The wall 
pressure ( pi ) is seen to change very rapidly at the point of the discontinuity. This large axial 
pressure gradient results in a large change in velocity at the discontinuity and thus the large change 
in skin friction coefficient that was seen in Fig. 3.5.7. 


For the smooth geometry the transverse pressure coefficient is shown in Fig. 3.5.9. The pressure 
changes are smooth throughout the entire length of the duct. The pressure coefficients shown in 
Fig. 3.5.9 are from calculations using a non-uniform grid. Fig. 3.5.10 shows identical results for 
calculations made using a uniform grid in the boundary layer. 


In Figs. 3.5.8, 3.5.9, and 3.5.10, a small transverse pressure gradient can be observed at the inlet to 
the duct. This small transverse pressure gradient is a result of specifying the inlet v-velocity to be 
equal to zero. This boundary condition would be typical of inviscid calculations ; however, the 
boundary layer in this viscous flow does have a finite v-velocity as a result of the growth of the 
boundary layer. The growth of the boundary layer for this test case is large because of the high 
absolute viscosity used. An additional calculation was made by modifying the boundary condition 
on the inlet v-velocity. From a previous calculation, the flow angles at the second axial station were 
recorded. These flow angles were then used as the inlet boundary condition for the inlet flow angle 
in a subsequent calculation. The transverse pressure coefficients for this case are shown in Fig. 
3.5.11. The inlet transverse pressure coefficient is reduced considerably when compared with the 
pressure coefficient shown in Fig. 3.5.9. When the Blasius solution was used to calculate the inlet 
v-velocity distribution, the solution over compensated and the transverse pressure gradient was in 
the opposite direction to that calculated when a zero v-velocity was specified. 
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It was also found to be important that the inlet and exit lengths of the straight portions of the duct 
were long enough. Oscillations in the static pressure would occur at the inlet and exit if these 
lengths were too short. 

An additional pressure coefficient was calculated for each of the cases described above. The pres- 
sure coefficient is defined as 


r = — 

0.5p« 1 2 


(3.5.2) 


where /jy,, is the freestream total pressure, p is the local static pressure, and 0.5pu 2 is the inlet dy- 
namic pressure. This pressure coefficient again shows the effect of curvature on the pressure 
through the boundary layer but it also shows the relative acceleration of the flow through the 
nozzle. The plots are 

1. non-smooth geometry with uniform grid ( Fig. 3.5.12 ) 


2. smooth geometry with non-uniform grid ( Fig. 3.5.13 ) 


3. smooth geometry with uniform grid (Fig. 3.5.14) 


4. inlet flow angle specified with non-uniform grid ( Fig. 3.5.15 ) 


stability and convergence 
the boundary layer. The 

1. [MV 1] no multi- volume used in a uniform grid 


This test case also provided us with an opportunity to investigate the 
properties of various multi- volume approaches (see section 2.8) within 
various strategies are 


3.0 RESULTS 


186 




A PI 

4* P3 
X P5 
O P7 
* P9 
X Pll 

VS X 


Fig. 3.5.12 


Pressure Coefficient 
Non-Smooth Geometry Uniform Grid 


3.0 RESULTS 


187 


PRESSURE COEFFICIENT 



A PI 

+ P3 
X P5 
<r> P7 
* P9 

VS X 


Fig. 3.5.13 Pressure Coefficient 

Smooth Geometry Non-Uniform Grid 


JS8 


3.0 RESULTS 



PRESSURE COEFFICIENT 



A PI 
+ P3 
X P 5 
O P7 
* P9 
X Pll 

VS X 


X/INLET HEIGHT 


Fig. 3.5.14 Pressure Coefficient 

Smooth Geometry Uniform Grid 


3.0 RESULTS 


189 



"1 



A pi 

+ P3 
X P5 
O P7 
+ P9 

VS X 


T 1 1 1 r 1 

0 12 3 4 5 

X/ I NILE T HEIGHT 

Fig. 3.5.15 Pressure Coefficient 

Non-Uniform Grid Inlet v-Velocity Specified 


3.0 RESULTS 


190 



2. [MV2] multi- volume used in a uniform grid with time steps based upon properties for the in- 
dividual control volume farthest away from the wall. 

3. [MV3] multi-volume used in a uniform grid with time steps based upon average properties for 
the entire multi-volume region. 

4. [MV4J multi-volume used in a non-uniform grid with time steps based upon properties for the 
individual control volume farthest away from the wall. 

5. [MV 5[ multi- volume used in a non-uniform grid with time steps based upon average properties 
for the entire multi-volume region. 

For the uniform and non-uniform grids, the multivolume approach was used for the 6 control 
volumes nearest to the wall. 

Two measures of convergence will be used in the following analysis. The maximum change in ve- 
locity in the flow field for one iteration divided by an average velocity for the flow field all multiplied 
by the time factor used in the time step determination will be our momentum residual, in other 
words, 

<u n + ' -u n ) 

momentum residual = S22L x TIMEF (3.5.3) 

u avg 

The time factor, TIMEF, is included so that calculations using different time factors in the time step 
determination can be judged on a common basis. A time factor of 4.0 was used for all of the laminar 
calculations. The mass flow rate error is defined as 

mass flow error = I ™ local — ^-| max x 100 (3.5.4) 

™in 

The momentum residual for multi-volume approaches 1 through 4 are shown in Fig. 3.5.16. The 
mass flow error for approaches 1 through 4 are shown in Fig. 3.5.17. The non-uniform grid shows 
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Fig. 3.5.16 Momentum Residual For Laminar Boundary 
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superior convergence properties when compared with the uniform grid results. When the uniform 
grid is used, approach MV3 is the best. 

When the momentum residual and mass flow error for the non-uniform grid approaches MV4 and 
MVS are compared in Figs. 3.5.18 and 3.5.19, the results are very similar. This agreement justifies 
the use of the simpler multi-volume approach ( MV4 ) for the non-uniform grid. The simpler 
multi-volume approach uses less computer time in the multi-volume part of the program. 

When a finer grid was used with the non-uniform grid arrangement in the boundary layer region 
(twice as many control volumes ), the skin friction coefficients were identical to those calculated 
from the coarse grid ( 5 points in the boundary layer). The momentum residual and the mass flow 
error for the coarse and fine grids are compared in Figs. 3.5.20 and 3.5.21 respectively. The con- 
vergence properties are very similar. 

The following conclusions can be reached from the above test case: 

1. the geometries should be inspected for discontinuities in the second derivative to ensure 
smooth solutions. 

2. uniform and non-uniform grids give essentially identical results for this laminar boundary layer. 

3. the non-uniform grid (with a factor of 2 change in spacing) has superior convergence proper- 
ties. 

4. the simple form of the multi-volume used with the non-uniform grid (MV4) is preferred be- 
cause of the reduced computational effort required. 
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DELTA V/ AVERAGE VELOCITY m TIMEF 
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Fig. 3.5.18 Momentum Residual For Laminar Boundary Layer 
Calculation Local Vs Global Time Steps 
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MASS FLOW ERROR 
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Fig. 3.5.19 Mass Flow Error For Laminar Boundary Layer 
Calculations Local Vs Global Time Steps 
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NUMBER OF ITERATIONS 


Fig. 3.5.20 Momentum Residual For Laminar Boundary Layer 
Calculations Coarse Vs Fine Grid 
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MASS FLOW ERROR 
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Fig. 3.5.21 Mass Flow Error For Laminar Boundary Layer 
Calculations Coarse Vs Fine Grid 
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3.6 SAMUEL AND JOUBERT INCOMPRESSIBLE 


TURBULENT BOUNDARY LA YER 


Incompressible turbulent boundary layer flow in a diverging duct was calculated for test case 0141 
of the Stanford Conference ( Samuel and Joubert ) (26). The geometry used by Moore (27) and 
the grid used in the present calculations are shown in Fig. 3.6.1. With this geometry, the top wall 
is treated as inviscid in the calculations. The inlet velocity is 26 m/s. The absolute viscosity is 
0.000018 kg/ m s. 

The duct height is 1 m at the inlet plane ( x = 0.85 m ). The length of the duct from the inlet to 
the exit plane is 3.3 m. The distribution of grid points across the duct and the inlet velocity profile 
are presented in Table 3.6.1. In the present calculations, the inlet total pressure is 102 kPa and the 
inlet total temperature is 300 K. The exit static pressure was chosen to be 100 kPa so that the inlet 
velocity in the freestream would be approximately 26 m/s. 

Fig. 3.6.2 shows a comparison of the calculated skin friction coefficient with the experimental re- 
sults and with the results from the Moore cascade flow program. The agreement is excellent. The 
skin friction coefficient, defined here as 


9 = 


o.5 P u; s . 


(3.6.1) 


uses a reference velocity , U„p of 27 m/s in its evaluation. 

A comparison of the calculated turbulent shear stress distribution , uV, with the experimental re- 
sults is shown is Figs. 3.6.3 and 3.6.4. The shear stresses are normalized with respect to the local 
freestream velocity. The agreement is very good. 


3.0 RESULTS 


199 




3.0 RESULTS 


200 


Table 3.6.1 Node Distribution and Inlet 

Velocity Profile for Samuel and Joubert 


1 

u/u 
e 

0.000165 

0,406 

0.000665 

0.574 

0.00200 

0.674 

0.00500 

0.760 

0.0110 

0.880 

0.0230 

0.990 

0.0470 

1.00 

0.0950 

- 

0.1910 

- 

0.3080 

- 

0.4145 

- 

0.5120 

- 

0.6275 

- 

0.7340 

- 

0.8405 

- 

0.9470 
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Fig. 3.6.3 Shear Stress Distributions at x = 1.04m, x = 1.44m, and 
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Pig . 3. 6. A Shear Stress Distributions at x = 2.38m, x = 2.89m, and 













Fig. 3.6.5 shows a comparison of the calculated and measured velocity profdes at two locations in 
the duct. The agreement is good at x= 2.87 m, however, the calculated boundary layer at x= 3.40 
m is thinner than the measurements. It is noteworthy that the program has been able to calculate 
an essentially incompressible flow ( M < 0.075 ). 


3.7 SAJ BEN'S DIFFUSER CALCULATIONS 


The diffuser geometry ( Model G ) (23) is shown in Fig. 3.7.1; the throat height ,h, was 44 mm and 
the ratio of the exit height to throat height was 1.513. The calculations begin at x/h= -3.6 and end 
at x/h = 8.2. Fig. 3.7. 1 also shows the computational grid used which had 87 grid points in the axial 
direction and 20 points across the flow. The radial distribution, of grid points used is shown in 
Table 3.7.1. The development of a turbulent boundary layer was modeled on both the curved and 
the flat walls. The inlet boundary layer thicknesses were specified as 9 % and 4.5 % of the inlet 
diffuser height for the curved and flat wall boundary layers, respectively. An absolute viscosity of 
0.000018 kg/ m s was used. The inlet total pressure was 135 kPa and the inlet total temperature 
was 300 K. 

For this calculation, the ratio of the exit static pressure to the inlet total pressure was 0.826. This 
computational exit static pressure is equal to the experimental exit static pressure plus a correction 
in pressure for the side wall boundary layer blockage. Suction slots upstream of the throat and 
downstream of the exit plane reduce the side wall boundary layer blockage and improve the two- 
dimensionality of the flow. However, the side wall boundary layers do affect the effective flow area 
of the diffuser and the two dimensional computations must reflect this when a suitable exit static 
pressure is chosen. Therefore the exit static pressure of 0.826 is 1.5 % greater than the experimental 
static pressure because of the blockage effect of the side wall boundary layers. In the experiment, 
this test point results in transonic flow in the diverging portion of the duct with a Mach number 
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Table 3.7.1 Transverse Distribution of Grid Points 


1 

0.0005 

2 

0.0020 

3 

0.0050 

4 

0.0110 

5 

0.0230 

6 

0.0470 

7 

0.0950 

8 

0.1910 

9 

0.3080 

10 

0.4305 

11 

0.5695 

12 

0.6920 

13 

0.8090 

14 

0.9050 

15 

0.9530 

16 

0.9770 

17 

0.9890 

18 

0.9950 

19 

0.9980 

20 

0.9995 
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of approximately 1.235 upstream of a nearly normal shock, and the flow remained fully-attached 
throughout the diffuser at this test condition. 

A plot of static pressure contours is shown in Fig. 3.7.2. The shock can be seen in the diverging 
portion of the duct. The shock is well defined as illustrated by the high clustering of contours at 
the shock. The M&M formula (see section 2.4) is used to calculate the the interpolated pressure 
which is used in the calculation of the effective density. Fig. 3.7.3 shows a Mach number contour 
plot for the calculations. The extent of the boundary layer can be seen from this figure. 

The calculated and measured curved wall static pressures are compared in Fig. 3.7.4. The shock is 
very well defined and no overshoot occurs in the static pressure. The static pressures do not agree 
so well downstream of the shock because the 2-D calculations do not reflect the rapid increase in 
the side wall boundary layer thickness and its effective blockage. The point of minimum static 
pressure in the calculations is located at xj h= 1.5. This is taken to be the location of shock. The 
Mach number upstream of the shock was determined to be 1.256 from the calculated total pressure 
ratio across the shock in the freestream. 

The computed and measured shock locations on the curved wall are compared in Fig. 3.7.5. The 
variable x ou is the shock location on the curved wall and x am is the shock location in the middle of 
the duct. The Mach number upstream of the shock, determined from static pressure measurements, 
is represented by in Fig. 3.7.5. 

The Mach number distribution through the nozzle at a fixed y/h of 0.0905 is shown in Fig. 3.7.6. 
This grid line was chosen because the maximum Mach number is located along it. The shock is 
sharp and no overshoots or undershoots occur. 

Comparisons of calculated and measured velocity profdes ( see Ref. 24) at four axial locations along 
the duct are shown in Figs. 3.7.7, 3.7.8, 3.7.9, and 3.7.10. The axial locations are x/h= 2.31, 4.03, 
6.34, and 8.2 respectively. The agreement is good especially at the two downstream stations. The 
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Fig. 3.7.4 Curved Wall Static Pressure For Sajben's Diffuser 
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Fig. 3.7.5 


Comparison of Computed and Measured Shock Position 
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A EXPERIMENT 
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Fig. 3.7.9 Velocity Profile at x/h = 6.34 
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velocities in Figs. 3.7.7-3.7.10 are normalized with respect to the maximum computed or measured 
velocity in the duct at that axial location, which ever is applicable. In the calculations, the edge of 
the boundary layer is located where the normalized total pressure gradient is 25.0. It may be noted 
that the present calculations give much better agreeement with the measured velocity profiles than 
the calculations of Liou, Coakley, and Bergmann of Reference 28 ( see Fig. 3.7.17). Other inves- 
tigators who have used Sajben's diffuser as a test case have not shown a comparison of their com- 
puted velocity profiles with the measurements of Ref. 24. 

Fig. 3.7.11 shows static pressure contours for calculations which were made using the three-point 
scheme for the interpolated pressure. The shock resolution is not nearly as well defined as that 
obtained when the M&M formula is used. Fig. 3.7.12 shows the corresponding Mach number 
contour plot when the three point interpolation scheme is used. 

The distribution of loss generation in the diffuser will be presented in three ways. First, the losses 
due to the curved wall boundary layer, the flat wall boundary layer, and the shock loss in the 
freestream are compared in Fig. 3.7.13. These losses were determined by first calculating the mass 
flow rates in the boundary layers and the freestream at the diffuser exit using the calculated 
boundary layer thicknesses as the boundaries between regions. Then the mass averaged total 
pressure loss at each axial location was calculated by integrating the total pressure loss out to these 
fixed mass flow rate values and then normalizing these losses with respect to the inlet freestream 
total pressure ( see Table 3.7.2). All three losses are approximately the same. The curved wall 
(bottom wall ) boundary layer contributes the largest proportion to the total losses because the 
boundary layer is thicker there. 

The total pressure loss along an inviscid streamline is compared with the mass averaged total pres- 
sure loss for the entire cross-section of the diffuser in Fig. 3.7. 14. This figure allows another means 
of comparing the shock and total losses. The total pressure loss through the shock is very well 
defined. The mass averaged total pressure loss through the shock is approximately 30% greater 
than the shock loss alone. 
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Table 3.7.2 Example Calculation of Total Pressure 
Loss for a Boundary Layer 


P . = ^ ^t-inlet ^t-local^ * 

t-losses — 

• 

m 

total 

* where it ^ is the total mass flow rate through 
total 

a given cross-section of the diffuser, and m j is the 
mass flow rate in the boundary layer at the exit. 
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The mass averaged total pressure divided by the inlet freestream total pressure at the diffuser exit 
is calculated from the numerical results to be 0.9615. From data given to us by M. Sajben and T. 
Bogar, the mass averaged total pressure calculated from the experimentally measured data is 0.965. 
The experimental data used for this calculation was measured midway between the side walls. The 
agreement is good. 

The measured maximum Mach number in the exit plane of 0.51 agrees well with the calculated 
value of .51 1. The boundary layer thicknesses are measured to be approximately 25% of the duct 
height on the curved wall and 23 % of the duct height on the flat wall. The calculations determined 
the boundary layer thickness on the curved wall to be 25% of the duct height and 23% of the duct 
height of the flat wall. 

The above calculations are after 5000 iterations using a TIMEF of 4.0. It may have been possible 
to run the calculation with a TIMEF of 2.0 after the initial transients but because of computer cost 
this option was not attemped for this test case. The momentum residual for this calculation is 
presented in Fig. 3.7.15. The continuity error after 5000 iterations was less than .1%. The unusual 
behavior observed between 3000-4000 iterations in Fig. 3.7.15 appears to be typical of the method 
in general when transonic flow in calculated. The same general convergence behavior is seen for 
the 1-D nozzle ( Ref. 12 ) in Fig. 3.7.16. A time factor of 2.0 was used for these 1-D calculations 
but the same behavior is seen here on a smaller time scale. The total CPU time for the Sajben 
calculations was approximately 35 minutes on the IBM 3031. 

A number of other workers have used this test case to verify the accuracy of their computational 
methods ( 28, 29, and 30 ). Liou, Coakley, and Bergmann's (28) calculations were made using a 
MacCormack type scheme in conjuction with a two-equation model for the turbulent stress pred- 
ictions. For the weak shock case, they used an exit static pressure of 0.8 * P lM „ . There was good 
agreement between their computed and measured wall static pressures. However, the calculated 
velocity profiles were not in good agreement with the experimental data. Fig. 3.7.17 shows a sample 
of their results. 
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Fig. 3.7.15 Momentum Residual For Sajben f s Diffuser 
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Fig. 3.7.17 Sample Results From Liou, Coakley, and Bergman 






Liu, Shamroth, and McDonald (29) used this test case to verify their computational method which 
solves the governing equations using a consistently split linearized block implicit scheme. A mixing 
length model was used to model turbulent stresses. For the weak shock case, an exit static pressure 
of 0.807 x P tJfxUt was used. They compared their computed wall static pressures with the exper- 
imental pressures and found good agreement (see Fig. 3.7.18). The variable ,a, referred to in Fig. 
3.7.18 is the value of the artificial dissipation parameter used in their calculations. They found that 
for values of a less than 0.1 the solution was invariant. Also presented in Ref. 28 were static 
pressure and Mach number contour plots for this weak shock case. These contour plots are similar 
to those in Figs. 3.7.2 and 3.7.3. The shock definition in the present work is sharper than Liu et. 
al. but fewer axial grid points were used in their calculations. No comparison was made by Liu, 
Shamroth, and McDonald between the measured and computed velocity profiles. 

Talcott and Kumar (30) also used the weak shock case for Sajben's diffuser as a test of their com- 
puter program's accuracy. They used an implicit MacCormack scheme to solve the governing 
equations and a Baldwin-Lomax mixing length turbulence model to predict the turbulent stresses. 
They also used an exit static pressure of 0.807 * P'M„r Fig. 3.7.19 shows a sample of their results. 
They present a static pressure contour plot for this weak shock case. As with Liu et. al., the shock 
definition of Talcott and Kumar's method is not as sharp as that obtained using the present method 
but again they have used fewer axial grid points in their calculations. Also, the static pressure 
contours at the throat of the diffuser are erratic. A plot of wall static pressure is also presented. 
No comparision was made between the computed and measured velocity profiles. 

In summary, the following observations can be made about the present calculations of the weak 
shock case in Sajben's diffuser. 

L There is better agreement between measured and calculated velocity profiles when the present 
method and static pressure specification are used. 
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Fig. 2 Pressure contours and velocity vector field in inlet for tur* 
bolent flow {P^/P l *0.933, M ( * 0.46). 



STSEAMWlSf LOCATION. »/* 'h = Mr<n> 

Fig. 3 Pressure distributions on the top surface of the inlet for tur- 
bulent flow (Af y * 0.46). 


Fig. 3.7.19 Sample Results From Talcott and Kumar 
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2. The present method compares favorably with other methods in the areas of static pressure and 
Mach number contour predictions and wall static pressure prediction. 

3. The agreement between calculated and measured mass averaged total pressure losses for the 
diffuser is excellent. 


3.8 ENERGY EQUA TION TEST CASES 


A number of test cases will be used to explore various aspects of the more complete form of the 
energy equation discussed in section 2.10 ( ENERGY EQUATION ). In order of presentation, the 
test cases are 


1. turbulent boundary layer flow in an adverse pressure gradient with an inlet ffeestream Mach 
number of 0.55. 


2. flat plate turbulent boundary layer with a freestream Mach number of 0.95. 

3. Sajben's diffuser calculations including the energy equation. 

4. flat plate turbulent boundary layer with a freestream Mach number of 2.8. 


Unless otherwise noted, the turbulent Prandtl number is a constant of 0.90 and the laminar Prandtl 


number is 0.73. The walls are adiabatic in all the calculations. 


TURBULENT BOUNDARY LAYER IN AN ADVERSE PRESSURE GRADIENT 

The geometry and grid used in this test case are identical to that used in Section 3.6 (Samuel and 
Joubert). The exit static pressure was reduced to 0.91 P t%in i, t so that the inlet freestream Mach 
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number was increased to 0.55. The same inlet velocity profile was specified as in section 3.6. The 
purpose of this test case was to illustrate the advantage of the new formulation of the energy 
equation that was discussed in Section 2.10. 


The work done due to pressure forces in the energy equation is 

work done due to pressure forces = — JJPu * • 

In the new formulation of the energy equation, the term 


+ irlfp u ' dA 


(3.8.1) 


was added to the calculation of the energy error for the transient solution. The purpose of including 
this additional term, Eq. 3.8.1, was to remove the continuity error contribution to the energy error 
caused by the work done due to pressure forces. If the base form of the energy equation 


p-^- x 5Ko/= - JJp ue t • dA + e,JJpy • dA - JJ - k VT- dA 

+ JJU# • (p Vu + \i~V~u T )\ ♦ dA - JJ/>« • dA . (3.8.2) 

is used in the current method, the transient solution results in a static temperature profile like that 
shown if Fig. 3.8.1 ( represented as + 's ). If the new form of the energy equation 

5v L 

p-£- x Wo1 = Y ~ JJp vbt'dA + A,JJpu • dA - JJ - kVT'dA 


+ JJ[“ • (F Vu + P Vu T )\ ♦ dA (3.8.3) 

is used, the transient solution results in a static temperature profile like that shown in Fig. 3.8.1 
(represented as triangles ). The results presented in Fig. 3.8.1 are from calculations after 500 iter- 
ations. It can clearly be seen that the new formulation gives a better transient solution to the energy 
equation and it should therefore result in a reduction in the computer time required to reach a 
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Fig. 3.8.1 Static Temperature Distribution Through The 

Boundary Layer For Samuel and Joubert Geometry 
at M=0.55, x=200 mm 
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steady state solution. Fig. 3.8.2 shows the corresponding total temperature profiles for the two 
formulations of the energy equation. 

FLAT PLATE BOUNDARY LAYER AT M = 0.95 


A turbulent boundary layer was calculated on the bottom wall of the geometry shown in Fig. 3.8.3. 
The height of the duct was 50 mm and the length of the duct was 200 mm. The computational 
grid shown in Fig. 3.8.3 consists of 21 axial grid points and 13 transverse grid points. The transverse 
distribution of grid points is shown in Table 3.8.1. The inlet boundary layer thickness was specified 
to be 10 % of the duct height. This grid placed 8 node points within the boundary layer. An ab- 
solute viscosity of 0.000018 kg/m s was used. The inlet total pressure was 105 kPa and the inlet 
total temperature was 300 K. The exit static pressure was 0.559 x P,M" which results in freestream 
Mach number of approximately 0.95. 

The total temperature distribution through the boundary layer at the exit of the duct is shown in 
Fig. 3.8.4. The total temperature at the near wall point ( y/h = 0.0001665 ) is 296.63 K which 
results in a recovery factor 


r = 


( Taw 7-qq) 

(Toco — Too) 


(3.8.4) 


of 0.927, where T„ is the adiabatic wall temperature, T„ is the freestream static temperature, and 
T 0ao is the freestream total temperature. The recovery factor for turbulent flow can be estimated 
by 

r = Pr x /3 . (3.8.5) 

For a Prandtl number of 0.73, the recovery factor is calculated from Eq. 3.8.5 to be 0.90. If an 
additional node point was located at y/h = 0.000055, the total temperature distribution remained 
essentially unchanged, however, the recovery factor changed to 0.910 because of the improved re- 
solution of the flow near the wall. To calculate the recovery factor with reasonable accuracy, it 
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Table 3.8.1 Transverse Distribution of Grid Points 
For M = 0.95 Test Case 


1 

0.0001665 

2 

0.000665 

3 

0.0020 

4 

0.005 

5 

0.011 

6 

0.023 

7 

0.047 

8 

0.09 5 

9 

0.191 

10 

0 . 308 

1 1 

0.4675 

12 

0.6805 

13 

0.8935 
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Fig. 3.8.4 Total Temperature Distribution For Flat Plate 
Boundary Layer at M = 0.95 
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was found that the near wall control volume boundary in the calculations needed to be at a distance 
ofj> + < 10. 


Beginning at the wall, the total temperature shown in Fig. 3.8.4 is first less them freestream total 
temperature and then in the outer part of the boundary layer the total temperature becomes greater 
than the freestream total temperature. Since energy should be conserved for this adiabatic wall case, 
the mass averaged difference between the freestream and local total temperature should be con- 
served in this boundary layer. The quantity, 


H = ^ P«(7q ~ Zo°o) 

PoO^OCjTqoo 


Ay 

~S~ 


(3.8.6) 


is plotted in Fig. 3.8.5. The residual does not go to zero at the freestream, however, the mass av- 
eraged total temperature for the boundary layer was .999987 T 0<o . 


When a laminar and turbulent Prandtl number of 1 was used in the calculations, the total temper- 
ature was essentially uniform throughout the boundary layer ( T„ = 0.99997T 0oo ). 


SAJBEN'S DIFFUSER CALCULATIONS INCLUDING THE ENERGY EQUATION 


Calculations were made in Section 3.7 of transonic flow in Sajben's diffuser with the assumption 
of constant total temperature. The calculations were restarted here at 5000 iterations with the en- 
ergy equation included. An additional 1 500 iterations were required for these calculations to con- 
verge sufficiently. The boundary conditions and grid are identical to those discussed in Section 3.7. 

A plot of static pressure contours is shown in Fig. 3.8.6. Fig. 3.8.7 shows a Mach number contour 
plot for the calculations. The results are essentially the same as those seen in Section 3.7 (see Figs. 
3.7.2 and 3.7.3). The velocity profiles at an x/h = 7.9 are compared in Fig. 3.8.8. The profiles are 
very similar and the assumption of constant total temperature results in a solution which is only 
slightly different than the calculations that were made using the full form of the energy equation. 
Total temperature profiles are presented in Fig. 3.8.9 at two axial locations ( x/h = 0.9 and x/h = 2.3) 
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Fig. 3.8.8 Sajben's Diffuser Velocity Profiles at x/h=7.9 
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Fig. 3.8.9 Total Temperature Distribution For Sajben's Diffuser 
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; in the diffuser. The freestream velocity is supersonic at x/h = 1.9 (before the shock) and subsonic 
at x/h = 3.5(afiter the shock). The calculated recovery factors of .95 and .98 are higher than the 
values of .93 and .97 expected from Eq. 3.8.5. The near wall control volume boundaries in these 
calculations were at y + = 32 and y + = 23 for the locations x/h = 1.9 and x/h = 3.2, respectively. 
These are greater than the minimum y + of 10 discussed earlier. 

FLAT PLATE TURBULENT BOUNDARY LAYER AT M = 2.8 

Van Driest (30) presents the total temperature distribution within a flat plate turbulent boundary 
layer with a freestream Mach number of 2.8. The experimental total temperature distribution is 
shown in Fig. 3.8.10. The geometry and grid for these calculations are shown in Fig. 3.8.1 1. The 
height of the duct was 63.5 mm and the length of the duct was 254 mm. The computational grid 
shown in Fig. 3.8.11 consists of 21 axial grid points and 14 transverse grid points. The transverse 
distribution of grid points is shown in Table 3.8.2. The inlet boundary layer thickness of 1/4 inch 
(6.35 mm) was 10 % of the duct height. The Reynolds number based upon axial distance is ap- 
proximately 10 7 . To stabilize these supersonic calculations, the upwind effective density method 
( see section 2.4 ) was used and the upstream boundary conditions were modified. The inlet velocity, 
total temperature, and total pressure were specified at the upstream boundary. Three calculations 
were performed with different assumptions about the turbulent Prandtl number. These assump- 
tions were 

1. Pr, = 0.90 Pr, = 0.73 

2. Pr, = 0.73 Pr, — 0.73 

3. Pr, varies linearly through the boundary layer from 0.9 at the wall to 0.66 in the freestream. 

The turbulent Prandtl number is typically set equal to a constant of 0.9 in calculations (21). The 
calculated total temperature distribution through this boundary layer using a constant turbulent 
Prandtl number of 0.9 is shown in Fig. 3.8.12 ( represented as triangles ). The recovery factor is 
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Fig. 24. Stagnation-temperature ratio variation acro>% the 
turbulent boundary layer at M x = 2.S. 


Fig. 3.8.10 Experimental Total Temperature Distribution 
for Van Driest 
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Table 3.8.2 Transverse Distribution of Grid Points 
For M = 2.8 Test Case 


1 

0.000055 

2 

0.00022 

3 

0.000665 

4 

0.0020 

5 

0.0050 

6 

0.0110 

7 

0.0230 

8 

0.0470 

9 

0.0950 

10 

0.1910 

11 

0.3080 

12 

0.4675 

13 

0.6805 

14 

0.8935 
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Fig. 3.8.12 Total Temperature Distribution For Flat Plate Boundary Layer at 


calculated to be 0.920 which compares with the empirically determined value of 0.90. However, the 
distribution of total temperature through the boundary layer does not compare well with the ex- 
periment. If the turbulent Prandtl number is set equal to the laminar Prandtl number of 0.73 , the 
total temperature distribution changes as seen in Fig. 3.8.12 ( represented as + 's ). The distribution 
through the outer part of the boundary layer has improved but the recovery factor of 0.813 does 
not compare well with the experimental value of 0.90. Schlichting (31) notes that the turbulent 
Prandtl number is not constant through the boundary layer. The experiments of H. Ludwieg (32) 
for turbulent flow through a pipe shows that the Prandtl number varies from approximately 0.9 at 
the pipe wall to 0.66 at the center of the pipe. This distribution is shown in Fig. 3.8.13. The var- 
iation is almost linear. For the third set of calculations, the Prandtl number was assumed to vary 
linearly through the boundary layer from 0.9 at the wall to 0.66 at the edge of the boundary layer. 
The total temperature distribution for this case is shown in Fig. 3.8.12 (represented as X's ). The 
total temperature distribution calculated using a variable Prandtl number is compared with the ex- 
perimental results in Fig. 3.8.14. Both the distribution of total temperature through the boundary 
layer and the recovery factor of 0.90 are in good agreement with the experimentally measured val- 
ues. 

The calculation of turbulent flow over a flat plate with a freestream Mach number of 0.95 was rerun 
using this variable turbulent Prandtl number and these results are compared in Fig. 3.8.15 with the 
total temperatures which were obtained with a constant Prandtl number of 0.9. The recovery factor 
was calculated to be 0.89 when a variable Prandtl number was used; this is in reasonable agreement 
with the value of 0.90 expected using Eq. 3.8.5. 
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Fig. 3.8.15 Total Temperature Distribution For Flat Plate 
Boundary Layer at M=0.95 With Variable Prandtl 
Number 
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4.0 SUMMARY AND CONCLUSIONS 


A finite volume time marching method has been extended to allow the calculation of laminar and 
turbulent flow in ducts. Both subsonic and supersonic flow can be calculated with the method. 
The starting point for the current method was the finite volume method of Denton ( 12 ). Denton's 
method is currently limited to inviscid flow. Viscous effects have been simulated by Denton and 
others using in viscid- viscous interaction approaches. The method that is presented in this disser- 
tation solves the Reynolds-averaged form of the Navier-Stokes equations directly. 

A number of features have been developed in the current work which allow the calculation of 
viscous flow with the finite volume time marching method. Some insight has also been shed upon 
the finite volume method in general in this work. 

FEATURES OF THE JINJTE VOLUME TIMELMARCHING METHOD 

The features which make the finite volume time marching method attractive can be summarized 
as 

1. The calculations are performed entirely in the physical domain. No coordinate transformations 
are needed. 
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2. Because of the finite volume formulation, mass and momentum are conserved for any grid. 


3. The explicit approach is attractive because of its simplicity. 

4. When the Euler or Navier Stokes equations are solved using the time marching approach, 
shocks develop naturally as part of the solution. 

ADVANCES MADE IN THE PRESENT WORK 

Many changes were needed to create a finite volume time marching method capable of calculating 

two dimensional laminar and turbulent flow. A much higher level of understanding of the method 

was needed to extend the finite volume time marching method into the viscous flow regime. Fea- 
tures which are new to this finite volume time marching method can be summarized as 

1. Control volumes were introduced which allow calculations to be made without smoothing of 
flow properties in the transverse direction. In particular, unlike Denton's method, there are the 
same number of control volumes as node points (same number of equations as unknowns). 

2. Convective forms of the momentum and energy equations were used to update the velocity 
and the total enthalpy. However, at the steady state, the governing equations returned to the 
conservative form. 

3. Time steps for use with the convective forms of the equations were derived. This resulted in 
time steps for the momentum and energy equations which are different from the time steps for 
the continuity equation. The time steps also vary spatially. 

4. The pressure is updated directly from the continuity equation rather than through the equation 
of state to aid in the calculation of pressure gradients in boundary layers. 
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5. A new pressure interpolation scheme was introduced which has improved shock capturing 
properties when compared with the standard 2 and 3 point interpolation schemes. Pressure 
interpolation is necessary for stability but here a stability criterion has been employed in de- 
riving a Mach number dependent interpolation scheme. 

6. Viscous forces that act on the control volume faces were calculated in the non-orthogonal 
physical grid arrangement. 

7. A multi-volume method for pressure changes in the boundary layer was introduced which al- 
lowed calculations to be made with very elongated control volumes. The multi-volume pre- 
serves the boundary layer approximation (static pressure changes through the boundary layer 
are small compared with the free stream velocity head). 

8. A method of discretization called transverse upwind differencing was developed to strengthen 
the centerpoint coefficient of the equations. This differencing enhances the stability of the 
method when the mass fluxes through the transverse faces become large compared with the 
mass fluxes through the stream wise faces. 

9. A new formulation for the energy equation was introduced which has improved transient be- 
havior when compared with the standard formulation. The new formulation removes the in- 
fluences of continuity and momentum errors from the energy equation during transients in the 
solution. 

JEST CASES, 

A number of test cases were used to illustrate the accuracy and the features of the new method. 

A large range of flow conditions were studied in these test cases . The freestream Mach numbers 

varied from 0.075 in the Samuel and Joubert test case ( section 3.6 ) to 2.8 in the Van Driest test 

case ( section 3.8 ). The results obtained in this dissertation can be summarized as follows. 
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1. FLAT PLATE LAMINAR BOUNDARY LAYER: The velocity profiles in a laminar flat 
plate boundary layer agreed well with the theoretical results. 

2. INVISCID CALCULATIONS OF SAJBEN'S DIFFUSER: The density and pressure up- 
date methods were compared here. When the density update method was used, the total 
pressure calculated using the effective pressure agreed with the theoretical distribution better 
than when the thermodynamic pressure was used. The total pressure distribution calculated 
using the effective density method agreed the best with the theoretical results. 

3. THE INFLUENCE OF TRANSVERSE SMOOTHING ON A STEP PROFILE IN A 
STRAIGHT DUCT: A step profile in inlet Mach number was drastically alterred by the 
transverse smoothing that is needed to stabilize the calculations using the old control volumes. 
When the new control volumes were used instead, the inlet Mach number profile changed very 
little through the duct. 

4. 1-D COMPUTATIONAL TESTS OF SHOCK CAPTURING USING PRESSURE IN- 
TERPOLATION FORMULAE TO CALCULATE EFFECTIVE DENSITY: Results from 
calculations using the 2 point, 3 point, and M&M interpolation formulas were compared with 
the theoretical 1-D solution through a converging-diverging nozzle. The new pressure in- 
terpolation scheme gave improved shock capturing properties when compared with the usual 
two and three point interpolation schemes. Overshoots and undershoots in properties were 
essentially eliminated when the Mach number dependent interpolation formula was used. All 
the interpolation formulae accurately predicted the overall total pressure loss through the 
shock. 

5. LAMINAR BOUNDARY LAYER IN TWO CONVERGING DUCTS: A laminar 

boundary layer was calculated on the curved wall of two converging ducts. One of ducts had 
a non-smooth geometry which resulted in solutions which were not smooth. In Fig 3.5.5, the 
skin friction coefficients from the finite volume calculations were compared with the skin fric- 
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tion coefficients calculated using the integral method of Thwaites and the agreement was good. 
A number of multivolume methods and grid arrangements were investigated and a non- 
uniform grid was found to have better convergence properties than a uniform grid. 

6. SAMUEL AND JOUBERT INCOMPRESSIBLE TURBULENT BOUNDARY LAYER: 
Incompressible turbulent boundary layer flow in an adverse pressure gradient was calculated 
using the current method. The agreement between the calculated and measured skin friction 
coefficient, turbulent shear stress distributions, and the mean velocity profiles was good. 

7. SAJBEN'S DIFFUSER CALCULATIONS: Transonic viscous flow through a converging 
diverging nozzle was calculated. The computed and measured velocity proflies were in good 
agreeement especially near the exit of the nozzle. The computed and measured shock location 
were compared in Fig. 3.7.6 and were found to be in good agreement. The ratio of PJP, <M „ 
was calculated to be .9615. This was in good agreement with the measured value of 0.965. 

8. ENERGY EQUATION TEST CASES: Four test cases were used to investigate introduction 
of the new formulation of the energy equation into the method. The new formulation was seen 
to give improved transient behavior when compared with the standard formulation. For tur- 
bulent flow over a flat plate with a freestream Mach number of 0.95, the calculated recovery 
factor of .91 compared well with the empirically determined value ( see Eq. 3.8.5 ) of 0.90. 
For flat plate boundary layer flow with a freestream Mach number of 2.8, the calculated total 
temperature profile was improved by using a variable Prandtl number through the boundary 
layer. The calculated recovery factor of 0.90 agreed very well with the empically determined 
value of 0.9. 

CONCLUSIONS 

The following important observations were made from these test cases and other experience with 

the current method: 
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1. The geometry of the duct used in the calculations must be smooth to obtain smooth solutions 
(see section 3.5). 

2. In the boundary layer region, a non-uniform grid in the transverse direction, with a factor of 
2 ratio in grid spacing, is preferred for stability reasons. A factor of 2 spacing is preferred be- 
cause the multi-volume procedure requires that pressure changes for each successively smaller 
multi-volume region decrease rapidly. 

3. If a factor of approximately 2 in grid spacing is used in the multi-volume region, the multi- 
volume procedure may be simplified by using the control volume furthest from the wall as the 
basis for the time step to be used for the entire multivolume region ( see section 3.5). 

4. The new pressure interpolation scheme gives improved shock definition when its results are 
compared with the results from the 2 and 3 point interpolation schemes (see section 3.4). 

5. The total pressure loss through a 1-D shock is predicted well by all of the interpolation 
schemes considered here(see section 3.4). 

6. The Prandtl mixing length model for the eddy viscosity gave good results for all of the viscous 
test cases investigated. 

7. The skin friction coefficients, velocity profiles, and turbulent shear stress distributions were 
calculated with good accuracy through the boundary layer with as few as 5 to 7 grid points in 
the boundary layer as was seen from the Samuel and Joubert results (section 3.6). 

8. Smoothing of flow properties, either linear or non-linear, as used by Denton, can greatly de- 
grade the accuracy of the results when large gradients in properties exist (see section 3.3). 

9. Transverse upwind differencing was needed to stabilize the transient solution of some of the 
test cases, however, there was no upwinding retained at the steady state. 
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10. The height of the control volume just outside of the multi- volume region should be approxi- 
mately equal to the sum of heights of the control volumes within the multivolume region. 
This arrangement is not necessary but is preferred so that pressure changes in the multi-volume 
and freestream regions match properly. 

11. A time factor of 4.0 is typically required during the initial part of a solution, but the time factor 
can typically be changed to 2.0 later in the solution. 

12. When the full energy equation was used, the adiabatic wall temperature was predicted well if 
a constant turbulent Prandtl number of 0.9 was used, however, the distribution of total tem- 
perature through the boundary layer was not predicted well (see section 3.8). 

13. The total temperature distribution in the outer part of the boundary layer was improved if a 
simple linear variation in Prandtl number was assumed through the boundary layer (see section 
3.8). 

14. The new formulation of the energy equation gives improved transient solution behavior (see 
section 3.8). 

15. The assumption of constant total temperature used in test cases 3. 1-3.7 should not adversely 
affect the accuracy of those results (see section 3.8). 

SUGGESTIONS FOR FURTHER WORK 

The current work is limited to attached flow. An important contibution would be to extend the 
method to be able to calculate regions of backflow. This would be especially important in the area 
of transonic flow since when the Mach number upstream of a normal shock exceeds approximately 
1.3 the flow will separate behind the shock. By the choice of control volumes, the flow direction 
in the present method must essentially be in the I grid direction. Therefore to permit the calculation 
of backflow, a different control volume arrangement may be needed. 
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The method could also be extended to calculate the flow through turbomachinery blade rows or 
blade rows as has been done by Dawes. The method may encounter problems at the leading and 
trailing edges of such geometries because of the long and thin control volume seen at the surfaces 
of the blades. By necessity, these long and thin control volumes must extend upstream and 
downstream of the leading and trailing edges, respectively. The current method has shown that it 
can calculate both laminar and turbulent flows, however, a model for transition would be needed 
if the method were to be applied to turbomachinery blade rows and cascades. The control volumes 
used in the present work have be extended to three dimensional inviscid flow (19), however, further 
work would be required to allow viscous calculations in three dimensions. Three dimensional ef- 
fects could be partially modelled by using a quasi-3-D approach; for example, the end wall 
boundary layer effects could be modelled by including an axial velocity density ratio in the calcu- 
lations. 

As can be seen, there is much more to investigate in this area but the present contribution has 
provided significant and unique contributions to the understanding of the time marching finite 
volume method in general and what is required to extend the method to viscous flow in particular. 
Hopefully, this work will provide insight to future workers in the finite volume area and will also 
inspire workers to try approaches which are different than the current computational philosophy. 
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Appendix A. TRUNCATION ERROR OF 
PRESSURE INTERPOLATION EQUATION 


The truncation error of the interpolated pressure used to calculate the density in Eq. 2.4.57 may 
be determined using a Taylor series analysis. The interpolated pressure P* is given by 

PU , = P, + ^(/> /+ ,-/>,) + a,(P /+ , - />,_ 0/2 + 02(P /+ , - P/_ 2 )/3 (Al.l) 

and to determine the accuracy of P* we will look at the magnitude of P‘, +l — P, +l . With grid 

spacing h, and expanding about 1+ 1, we have 

P I - 2 = P ~ 3 hP' + 9(A 2 /2)p" - 0(A 3 ) (/1.1.2a) 

Pi - , = P - 7 hP' + A(h 2 j2)p" - 0(h 3 ) (A A. 2b) 

P,= P- hF + (h 2 l2)p" - 0(h 3 ) (/1.1.2c) 

P/ +1 = P (AA.2d) 

Therefore, 
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P e i+\ -/»/+,= Koq + a t + 02 - \)P' - (h 2 /2)(ao + 2a x + 3^ - 1 )p" + 0(h 3 ). 04.1.3) 


And if 


Oq + a, + aj = 1 


04 . 1 . 4 ) 


then the difference between P* and P is of the order of h 1 , so that P' is a second order accurate 
approximation for P. 
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Appendix B. VELOCITY GRADIENT IN A 
NON-ORTHOGONAL GRID 


To calculate the velocity gradient, Vu Lt for a nonorthogonal grid system, we must be able to 
transform gradients in the grid directions (14, and K) into gradients in the coordinate directions 
(x,y, and z). From the transformation law of a vector, the gradient of a scalar in one coordinate 
system can be related to the gradient in another coordinate system using the direction cosines be- 
tween the two coordinate systems. We will use the cartesian system as our base coordinate system 
here. Let and £) K be directional vectors which have a direction parallel with the non- 
orthogonal coordinate system and span a unit change in 14, or K respectively. These vectors are 
shown in Fig. A 2.1. They can be described in terms of their cartesian components. In other 
words, we can represent them as 


£/ = a x i + btf + c,& 

{A. 2.1) 

Dj = <*21 + bil+ c 2 k. 

{A. 2.2) 

D. k = + bjl + c 3 & 

(A. 2.3) 
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where j,^’, and £ are unit vectors in the cartesian coordinate directions, and a,b, and c are the 
magnitudes of the components of the directional vectors. 

The change of a variable ,u, from one point to another within a grid system can be represented as 

-f ^ = u /+1 - W/ = J, /+1 cfc-Fu 04.2.4) 

where r is a position vector between the two nodes. Rather than dealing with a local derivative in 
a continuum , we are dealing with approximations of derivatives using finite changes. If the direc- 
tional vector, £>/, is chosen such that its termini correspond to the nodes at which the properties 
are evaluated in Eq. A.2.4, and taking Vu as uniform between grid points, 

= (A. 2.5) 


But D.i = £>+ 1 - Hi therefore, 


du 

dl 


= Hr Vu 


(A. 2.6) 



Vu 


(A. 2.7) 


lh c ‘' F “ 


(A, 2.8) 


where IJ, and K are the coordinate indices for the non-orthogonal coordinate system. Expanding 
equations A 2.6, A 2.7, and A 2.8, we get 


du _ 
dl 


a \ 


du * l du 
dx oy 


+ 



(A. 2.9) 


= cJlL + bto + c.i“ 
3J ^dx 02 dy 2 dz 


(/4.2.10) 
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(-4.2-11) 



fi- + *,4^ + 

ox oy 



where a, = D u , the x component of D,\ £>, — D Tj etc... Equations A.2.7, A.2.10, and A.2.I1 can 
be combined into a single matrix equation, 


r du\ 

<*i 

bx 

c l 


{du ' 

dl 


dx 

du ( _ 


i 



du 

1 dJ l 


h 


* 

by 

| du \ 


L 



du 

dKJ 

_ a 3 

b 3 

c 3 

1 

\dzj 


(A. 2.12) 


For the calculation of viscous stresses, Eq. 2.7.5, we need the gradients of velocity ,Vu L , in cartesian 
coordinates. To obtain this gradient, we need to take the inverse of the above matrix [A], in other 
words, 




04.2.13) 


The inverse of matrix A , [A ]~ ', will now be determined. From linear algebra, 


[A\~ l = -±-{adj\A}) 


(A. 2.14) 


where 1/4 1 is the determinant of (/4) and adj[/4] is the adjoint of [A]. The determinant of [A] is 


det[/4] = 


a, 6, c, 
<h h c i 


a 3 h c 3 


a l 


b 2 

b 3 c 3 


~ 


°2 <>1 


a 3 c 3 


+ C, 


°2 *2 
a 3 *3 
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a l(*2 c 3 ~ h 0 !) ~ ^l(«2 c 3 “ a 3C2) + Cl(<*2*3 ~ «A) • 
This can be identified as 

det[/t] = D[ • (Dj x D^) 

Now let us represent the adjoint of A , adj(,4], in the form 

^ll ^12 ^13 
adj[A\ = /l 2 i ^22 ^23 
^31 ^32 ^33 

where the components of the adjoint matrix are 

bj c 2 

^11 “ ~ ® 2 C 3 ~ 

*3 c 3 


(,4.2.15) 


(>4.2.16) 


(>4. 2. 17) 


(A. 2. 18) 


However, >4 n , can be identified as the x-component of Qj x D K , in other words, 

^2 C 3 — ^3 C 2 ~ x component of Dj x J) K (,4.2.19) 

The components of the adjoint matrix, A 12 and A l3 , are the y and z components of Di x 5*. Both 
the results from the numerator and denominator of Eq. A.2.12, lead to the result that 


Dj x Dx + Dk x Di + D[ x Dj Su L 
Z2[ ’ (Dj x Dk) dl 22/ • (22/ x 22^) 57 22/ • (22j x 22^-) 


(,4.2.20) 


This result can be generalized to the other two components of velocity, v and w. This is the result 
that is used to calculate velocity gradients in the current program. 
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